Coverage for python/src/dolfinx_mpc/multipointconstraint.py: 91%
372 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-07 09:58 +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, Sequence, 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
17import ufl
19import dolfinx_mpc.cpp
21from .container import (
22 _UNSET,
23 MPCData,
24 _cpp_function,
25 _deprecated,
26 _float_array_types,
27 _float_classes,
28 _mpc_classes,
29 _mpc_data_classes,
30 _scalar_type,
31 _tolerance,
32 _Unset,
33)
34from .dictcondition import create_dictionary_constraint
35from .integralcondition import create_integral_constraint
36from .rbe import create_rbe2, create_rbe3
39class MultiPointConstraint:
40 """
41 Hold data for multi point constraint relation ships,
42 including new index maps for local assembly of matrices and vectors.
44 The constraint is affine, :math:`x = K x_{red} + g`, where :math:`g` is
45 supplied through `rhs_coeffs` and through the Dirichlet conditions in
46 `bcs`. With neither, :math:`g=0` and the constraint is the usual linear
47 one.
49 Args:
50 V: The function space
51 dtype: The scalar type of the coefficients, used by everything built from this
52 constraint. Defaults to the default scalar type of DOLFINx, real or complex, at the
53 precision of the mesh of `V`. Its precision must be that of the mesh.
54 bcs: Dirichlet boundary conditions for the problem. A master degree of
55 freedom that is constrained by one of these is removed from the
56 equation of its slave, and its contribution folded into the
57 constraint offset :math:`g`. As the offset is recomputed from the
58 current values of the conditions by :func:`update_constants`, time
59 dependent boundary data is supported.
60 rhs_coeffs: Function holding an additional inhomogeneity :math:`g_s`
61 for the slave degrees of freedom, i.e.
62 :math:`u_s = \\sum_j c_j u_{m_j} + g_s`.
63 """
65 _slaves: npt.NDArray[numpy.int32]
66 _masters: npt.NDArray[numpy.int64]
67 _coeffs: _float_array_types
68 _owners: npt.NDArray[numpy.int32]
69 _offsets: npt.NDArray[numpy.int32]
70 _master_spaces: List[tuple[npt.NDArray[numpy.int32], Optional[_fem.FunctionSpace]]]
71 _bcs: List[_fem.DirichletBC]
72 _rhs_coeffs: Optional[_fem.Function]
73 _scale_function: Optional[_fem.Function]
74 _rbe2: List[list]
75 _rbe3: List[tuple]
76 _rbe3_data: Optional[tuple]
77 V: _fem.FunctionSpace
78 _input_space: _fem.FunctionSpace
79 finalized: bool
80 _cpp_object: _mpc_classes
81 _dtype: npt.DTypeLike
82 __slots__ = tuple(__annotations__)
84 def __init__(
85 self,
86 V: _fem.FunctionSpace,
87 dtype: npt.DTypeLike | None = None,
88 bcs: Optional[List[_fem.DirichletBC]] = None,
89 rhs_coeffs: Optional[_fem.Function] = None,
90 ):
91 dtype = _scalar_type(V.mesh.geometry.x.dtype, dtype)
92 self._slaves = numpy.array([], dtype=numpy.int32)
93 self._masters = numpy.array([], dtype=numpy.int64)
94 self._coeffs = numpy.array([], dtype=dtype) # type: ignore
95 self._owners = numpy.array([], dtype=numpy.int32)
96 self._offsets = numpy.array([0], dtype=numpy.int32)
97 self._master_spaces = []
98 self._bcs = [] if bcs is None else list(bcs)
99 if rhs_coeffs is not None:
100 if not rhs_coeffs.x.array.dtype == dtype:
101 raise ValueError("rhs_coeffs must have the same dtype as the MPC")
102 if rhs_coeffs.function_space != V:
103 raise ValueError("rhs_coeffs must be a Function in the space of the constraint")
104 self._rhs_coeffs = rhs_coeffs
105 self._scale_function = None
106 # Per space on a spider mesh: [W, the tied space, the block of W], the block set by finalize
107 self._rbe2 = []
108 # Feet of RBE3 constraints on this space, (V, dofs, spiders, weights) per call, built by
109 # finalize, which keeps the arrays for update_rbe3
110 self._rbe3 = []
111 self._rbe3_data = None
112 self.V = V
113 # Kept after finalize replaces `V` by the extended space, which contains no Dirichlet condition
114 self._input_space = 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 master_space: Optional[_fem.FunctionSpace] = None,
127 master_blocks: Optional[npt.NDArray[numpy.int32]] = None,
128 ):
129 """
130 Add new constraint given by numpy arrays.
132 Args:
133 V: The function space for the constraint
134 slaves: List of all slave dofs (using local dof numbering) on this process
135 masters: List of all master dofs (using global dof numbering) on this process
136 coeffs: The coefficients corresponding to each master.
137 owners: The process each master is owned by.
138 offsets: Array indicating the location in the masters array for the i-th slave
139 in the slaves arrays, i.e.
141 .. highlight:: python
142 .. code-block:: python
144 masters_of_owned_slave[i] = masters[offsets[i]:offsets[i+1]]
146 master_space: The function space all masters belong to, if not `V`. It must be the
147 space of another constraint finalized together with this one by
148 :func:`finalize_multipointconstraints`, or a subspace of it, and `masters` is in
149 the global numbering of that constraint's space.
150 The masters of a slave may then be in another block of a blocked problem.
151 master_blocks: The block of each master, for masters from several spaces: its position
152 in the list of constraints given to :func:`finalize_multipointconstraints`. Each
153 master is in the global numbering of its block. Exclusive with `master_space`.
155 Note:
156 Collective when `master_space` or `master_blocks` is given: every process must call
157 it with the same `master_space`, or with `master_blocks` (possibly empty).
158 """
159 assert V == self.V
160 self._raise_if_finalized()
161 if master_space is not None and master_blocks is not None:
162 raise ValueError("Give either master_space or master_blocks, not both")
163 if master_blocks is not None and len(master_blocks) != len(masters):
164 raise ValueError("master_blocks must have one entry per master")
166 # Recorded on every process, also without local slaves, so that every process resolves the
167 # blocks of the masters the same way when the constraints are finalized
168 if master_blocks is not None:
169 self._master_spaces.append((numpy.asarray(master_blocks, dtype=numpy.int32), None))
170 elif master_space is not None:
171 self._master_spaces.append((numpy.full(len(masters), -1, dtype=numpy.int32), master_space))
172 else:
173 self._master_spaces.append((numpy.full(len(masters), -2, dtype=numpy.int32), None))
174 if len(slaves) > 0:
175 self._offsets = numpy.append(self._offsets, offsets[1:] + len(self._masters))
176 self._slaves = numpy.append(self._slaves, slaves)
177 self._masters = numpy.append(self._masters, masters)
178 self._coeffs = numpy.array(numpy.append(self._coeffs, coeffs), dtype=self._dtype)
179 self._owners = numpy.append(self._owners, owners)
181 def add_integral_constraint(
182 self,
183 weight_form,
184 value,
185 bcs: Optional[List[_fem.DirichletBC]] = None,
186 rtol: numpy.floating | float | None = None,
187 *,
188 coefficient_tol: Optional[float] = None,
189 ):
190 r"""Constrain a scalar integral of the solution, :math:`L(u) = \gamma`.
192 The functional is given as a linear form, and turned into a constraint
193 with a single slave by :func:`dolfinx_mpc.create_integral_constraint`;
194 see there for the derivation and the cost. The inhomogeneity
195 :math:`\gamma/w_s` is written into the ``rhs_coeffs`` function of this
196 constraint, which is created here if none was supplied to the
197 constructor.
199 Args:
200 weight_form: A linear form in ``ufl.TestFunction(V)`` defining the
201 functional, for instance ``v * ufl.dx``. Its test function must
202 be in the function space of this constraint.
203 value: The prescribed value :math:`\gamma` of the functional.
204 bcs: Dirichlet conditions on the space. A constrained degree of
205 freedom is never chosen as the slave. Pass the same conditions
206 to the constructor to have a constrained *master* folded into
207 the constraint offset. Defaults to the conditions given to the
208 constructor.
209 rtol: Deprecated, use `coefficient_tol`.
210 coefficient_tol: A master whose coefficient is below
211 `coefficient_tol` times the largest one is dropped. Defaults to
212 `500` machine epsilon of the real type of the constraint.
214 Note:
215 Collective. Must be called by every process.
216 """
217 self._raise_if_finalized()
218 if rtol is not None:
219 _deprecated("rtol", "`coefficient_tol`")
220 coefficient_tol = rtol if coefficient_tol is None else coefficient_tol
221 slaves, masters, coeffs, owners, offsets, rhs = create_integral_constraint(
222 self.V,
223 weight_form,
224 value,
225 self._bcs if bcs is None else bcs,
226 coefficient_tol=_tolerance(coefficient_tol, self._dtype),
227 )
228 if self._rhs_coeffs is None:
229 self._rhs_coeffs = rhs
230 else:
231 # Slaves of separate constraints are disjoint, so the offsets add
232 self._rhs_coeffs.x.array[:] += rhs.x.array
233 self.add_constraint(self.V, slaves, masters, coeffs, owners, offsets)
235 def add_constraint_from_mpc_data(
236 self,
237 V: _fem.FunctionSpace,
238 mpc_data: Union[_mpc_data_classes, MPCData],
239 master_space: Optional[_fem.FunctionSpace] = None,
240 ):
241 """
242 Add new constraint given by an `dolfinc_mpc.cpp.mpc.mpc_data`-object. See
243 :meth:`add_constraint` for `master_space`.
244 """
245 self._raise_if_finalized()
246 self.add_constraint(
247 V,
248 mpc_data.slaves,
249 mpc_data.masters,
250 mpc_data.coeffs,
251 mpc_data.owners,
252 mpc_data.offsets,
253 master_space=master_space,
254 )
256 def finalize(self, filter: Optional[numpy.floating] = None) -> None:
257 """
258 Finializes the multi point constraint. After this function is called, no new constraints can be added
259 to the constraint. This function creates a map from the cells (local to index) to the slave degrees of
260 freedom and builds a new index map and function space where unghosted master dofs are added as ghosts.
262 Args:
263 filter: If given, discard every master whose coefficient satisfies
264 :math:`|c_{sj}| < \\mathrm{filter}\\cdot\\max_k|c_{sk}|`, the
265 maximum being over the masters of that same slave. A negligible
266 coefficient contributes nothing to the constraint, but still
267 costs a ghost, a row of the sparsity pattern and an entry in
268 every element matrix modification, so removing them can shrink
269 :math:`K^HAK` substantially. With `None` (the default) every
270 master supplied is kept.
272 Note:
273 Filtering changes the constraint that is enforced, by exactly the
274 terms that are dropped. It is local and adds no communication.
276 Note:
277 To finalize the constraints of several function spaces, for instance the blocks of a
278 :class:`ufl.MixedFunctionSpace`, use :func:`finalize_multipointconstraints`.
279 """
280 finalize_multipointconstraints([self], filter)
282 def update_constants(self) -> None:
283 """
284 Recompute the constraint offset :math:`g` from the current values of the Dirichlet
285 conditions supplied to the constructor.
287 Call this whenever the value of one of those conditions changes, for instance between
288 time steps, before re-assembling. :class:`LinearProblem` calls it automatically.
290 Note:
291 Collective. Must be called by every process.
292 """
293 self._raise_if_not_finalized()
294 if self._rhs_coeffs is not None:
295 # Pass the array natively. Zero-copy, zero-allocation.
296 num_dofs_local = self.V.dofmap.index_map_bs * (
297 self.V.dofmap.index_map.size_local + self.V.dofmap.index_map.num_ghosts
298 )
299 rhs_coeffs = self._rhs_coeffs.x.array[:num_dofs_local]
300 self._cpp_object.set_rhs_coeffs(rhs_coeffs)
302 self._cpp_object.update_constants()
304 @property
305 def constants(self) -> _float_array_types:
306 """
307 The constraint offset :math:`g` for each degree of freedom local to the process,
308 i.e. the affine term in :math:`x = K x_{red} + g`.
309 """
310 self._raise_if_not_finalized()
311 return self._cpp_object.constants
313 @property
314 def has_inhomogeneity(self) -> bool:
315 """
316 Whether any process carries a non-zero constraint offset. The value is globally
317 reduced, so it is identical on every process.
318 """
319 self._raise_if_not_finalized()
320 return self._cpp_object.has_inhomogeneity
322 @property
323 def master_blocks(self) -> npt.NDArray[numpy.int32]:
324 """
325 The block of each master, parallel to ``masters.array``: the position, in the list given
326 to :func:`finalize_multipointconstraints`, of the constraint whose space the master is in.
327 The local index of a master is in the space of its block.
328 """
329 self._raise_if_not_finalized()
330 return self._cpp_object.master_blocks
332 @property
333 def dtype(self) -> type:
334 """The scalar type of the coefficients, also that of everything built from the constraint."""
335 return self._dtype
337 @property
338 def has_cross_block_masters(self) -> bool:
339 """Whether a master on any process is in another block than the slaves."""
340 self._raise_if_not_finalized()
341 return self._cpp_object.has_cross_block_masters
343 def create_periodic_constraint_topological(
344 self,
345 V: _fem.FunctionSpace,
346 meshtag: _mesh.MeshTags,
347 tag: int,
348 relation: Callable[[numpy.ndarray], numpy.ndarray],
349 bcs: List[_fem.DirichletBC],
350 scale: Union[_float_classes, float, complex] = 1.0,
351 tol: Union[_float_classes, float, None, _Unset] = _UNSET,
352 num_threads: Optional[int] = 1,
353 *,
354 distance_tol: Optional[float] = None,
355 coefficient_tol: Optional[float] = None,
356 ):
357 """
358 Create periodic condition for all closure dofs of on all entities in `meshtag` with value `tag`.
359 :math:`u(x_i) = scale * u(relation(x_i))` for all of :math:`x_i` on marked entities.
361 Args:
362 V: The function space to assign the condition to. Should either be the space of the MPC or a sub space.
363 meshtag: MeshTag for entity to apply the periodic condition on
364 tag: Tag indicating which entities should be slaves
365 relation: Lambda-function describing the geometrical relation
366 bcs: Dirichlet boundary conditions for the problem (Periodic constraints will be ignored for these dofs)
367 scale: Factor of the masters, of the scalar type of the constraint
368 tol: Deprecated, use `distance_tol` and `coefficient_tol`: a value sets both, `None` sets
369 `coefficient_tol=0`.
370 num_threads: The number of threads to use for certain operations
371 distance_tol: The largest distance from a mapped slave point to a master cell for the point to
372 be in the cell, and the padding of the bounding boxes of the cells. Defaults to `500`
373 machine epsilon of the coordinate type of the mesh.
374 coefficient_tol: A master whose coefficient is below `coefficient_tol` times the largest of
375 its slave is dropped. `0` keeps every master, so that the coefficients can later be changed
376 with :func:`scale_coefficients` or :func:`update_coefficients`. Defaults to `500`
377 machine epsilon of the real type of the constraint.
378 """
379 bcs_ = [bc._cpp_object for bc in bcs]
380 if isinstance(scale, numpy.generic): # nanobind conversion of numpy dtypes to general Python types
381 scale = scale.item() # type: ignore
382 distance_tol, coefficient_tol = self._tolerances(distance_tol, coefficient_tol, tol=tol)
383 is_input_space = V is self.V
384 if not (is_input_space or self.V.contains(V)):
385 raise RuntimeError("The input space has to be a sub space (or the full space) of the MPC")
386 mpc_data = _cpp_function("create_periodic_constraint_topological", self._dtype)(
387 V._cpp_object,
388 meshtag._cpp_object,
389 tag,
390 relation,
391 bcs_,
392 scale,
393 not is_input_space,
394 distance_tol,
395 coefficient_tol,
396 num_threads,
397 )
398 self.add_constraint_from_mpc_data(self.V, mpc_data=mpc_data)
400 def create_periodic_constraint_geometrical(
401 self,
402 V: _fem.FunctionSpace,
403 indicator: Callable[[numpy.ndarray], numpy.ndarray],
404 relation: Callable[[numpy.ndarray], numpy.ndarray],
405 bcs: List[_fem.DirichletBC],
406 scale: Union[_float_classes, float, complex] = 1.0,
407 tol: Union[_float_classes, float, None, _Unset] = _UNSET,
408 num_threads: Optional[int] = 1,
409 *,
410 distance_tol: Optional[float] = None,
411 coefficient_tol: Optional[float] = None,
412 ):
413 """
414 Create a periodic condition for all degrees of freedom whose physical location satisfies
415 :math:`indicator(x_i)==True`, i.e.
416 :math:`u(x_i) = scale * u(relation(x_i))` for all :math:`x_i`
418 Args:
419 V: The function space to assign the condition to. Should either be the space of the MPC or a sub space.
420 indicator: Lambda-function to locate degrees of freedom that should be slaves
421 relation: Lambda-function describing the geometrical relation to master dofs
422 bcs: Dirichlet boundary conditions for the problem
423 (Periodic constraints will be ignored for these dofs)
424 scale: Factor of the masters, of the scalar type of the constraint
425 tol: Deprecated, use `distance_tol` and `coefficient_tol`: a value sets both, `None` sets
426 `coefficient_tol=0`.
427 num_threads: The number of threads to use for certain operations.
428 distance_tol: The largest distance from a mapped slave point to a master cell for the point to
429 be in the cell, and the padding of the bounding boxes of the cells. Defaults to `500`
430 machine epsilon of the coordinate type of the mesh.
431 coefficient_tol: A master whose coefficient is below `coefficient_tol` times the largest of
432 its slave is dropped. `0` keeps every master, so that the coefficients can later be changed
433 with :func:`scale_coefficients` or :func:`update_coefficients`. Defaults to `500`
434 machine epsilon of the real type of the constraint.
435 """
436 if isinstance(scale, numpy.generic): # nanobind conversion of numpy dtypes to general Python types
437 scale = scale.item() # type: ignore
438 distance_tol, coefficient_tol = self._tolerances(distance_tol, coefficient_tol, tol=tol)
439 bcs = [] if bcs is None else [bc._cpp_object for bc in bcs]
440 is_input_space = V is self.V
441 if not (is_input_space or self.V.contains(V)):
442 raise RuntimeError("The input space has to be a sub space (or the full space) of the MPC")
443 mpc_data = _cpp_function("create_periodic_constraint_geometrical", self._dtype)(
444 V._cpp_object,
445 indicator,
446 relation,
447 bcs,
448 scale,
449 not is_input_space,
450 distance_tol,
451 coefficient_tol,
452 num_threads,
453 )
454 self.add_constraint_from_mpc_data(self.V, mpc_data=mpc_data)
456 def create_submesh_constraint(
457 self,
458 V: _fem.FunctionSpace,
459 master_space: _fem.FunctionSpace,
460 entity_map: _mesh.EntityMap,
461 bcs: Optional[List[_fem.DirichletBC]] = None,
462 scale: Union[_float_classes, float, complex] = 1.0,
463 tol: Union[_float_classes, float, None, _Unset] = _UNSET,
464 num_threads: int = 1,
465 *,
466 coefficient_tol: Optional[float] = None,
467 ):
468 r"""
469 Tie the degrees of freedom of `V` to `master_space` on a related mesh: a submesh and its
470 parent, related by `entity_map` as returned by :func:`dolfinx.mesh.create_submesh`.
472 Every degree of freedom of `V` in the closure of a cell related to a cell of
473 `master_space` becomes a slave, :math:`u(x_i) = \mathrm{scale}\, u_m(x_i)`, with
474 :math:`u_m` evaluated in the related cell and component `b` tied to component `b`. With `V`
475 on the submesh this is every degree of freedom of `V`, for instance the trace
476 :math:`\bar u = u|_\Gamma` on a submesh of facets; with `V` on the parent it is the
477 degrees of freedom on the submesh. No search is involved: the related cell is a table
478 lookup. For a submesh of facets the parent cell is one attached to the facet, so for a
479 discontinuous `master_space` the side is arbitrary.
481 Args:
482 V: The space of the constraint, or a subspace of it
483 master_space: The space of another constraint finalized together with this one by
484 :func:`finalize_multipointconstraints`, or a subspace of it. Its mesh and the mesh
485 of `V` are the two meshes of `entity_map`, either way round.
486 entity_map: Relates the cells of the submesh to entities of the parent, of
487 codimension 0 or 1
488 bcs: Dirichlet conditions on the space of the constraint. Their degrees of freedom
489 are not made slaves.
490 scale: Factor of the masters, of the scalar type of the constraint
491 tol: Deprecated, use `coefficient_tol`: a value sets it, `None` sets `coefficient_tol=0`.
492 num_threads: The number of threads to use
493 coefficient_tol: A master whose coefficient is below `coefficient_tol` times the largest of
494 its slave is dropped. `0` keeps every basis function of the related cell, so that the
495 coefficients can later be changed with :func:`scale_coefficients` or
496 :func:`update_coefficients`. Defaults to `500` machine epsilon of the real type of the
497 constraint. No distance tolerance is needed, as the related cell is not searched for.
499 Raises:
500 ValueError: If `entity_map` does not relate the two meshes, relates entities other
501 than the cells of the submesh, is of codimension above 1, or the spaces have
502 different numbers of components. Raised on every process.
504 Note:
505 Collective.
506 """
507 self._raise_if_finalized()
508 if not (V is self.V or self.V.contains(V)):
509 raise ValueError("V must be the space of the constraint or a subspace of it")
510 if isinstance(scale, numpy.generic): # nanobind conversion of numpy dtypes to general Python types
511 scale = scale.item() # type: ignore
512 if not isinstance(tol, _Unset):
513 _deprecated("tol", "`coefficient_tol`")
514 if coefficient_tol is None:
515 coefficient_tol = 0.0 if tol is None else tol
516 bcs_ = [] if bcs is None else [bc._cpp_object for bc in bcs]
517 mpc_data = _cpp_function("create_submesh_constraint", self._dtype)(
518 V._cpp_object,
519 master_space._cpp_object,
520 entity_map._cpp_object,
521 bcs_,
522 scale,
523 _tolerance(coefficient_tol, self._dtype),
524 num_threads,
525 )
526 self.add_constraint_from_mpc_data(self.V, mpc_data=mpc_data, master_space=master_space)
528 def _add_rbe2(self, dofs: list[npt.NDArray[numpy.int32]], W: _fem.FunctionSpace, x=None):
529 """Tie `dofs[k]` to spider `k` of `W`, and record `W` for :meth:`update_rbe2`."""
530 spiders = [numpy.full(len(d), k, dtype=numpy.int64) for k, d in enumerate(dofs)]
531 mpc_data = create_rbe2(
532 self.V,
533 numpy.concatenate(dofs) if dofs else numpy.zeros(0, dtype=numpy.int32),
534 numpy.concatenate(spiders) if spiders else numpy.zeros(0, dtype=numpy.int64),
535 W,
536 self._dtype,
537 x,
538 )
539 self.add_constraint_from_mpc_data(self.V, mpc_data=mpc_data, master_space=W)
540 if not any(W is entry[0] for entry in self._rbe2):
541 self._rbe2.append([W, self.V, None])
543 def add_rbe2_topological(
544 self,
545 dim: int,
546 entities: Union[npt.NDArray[numpy.int32], Sequence[Optional[npt.NDArray[numpy.int32]]]],
547 W: _fem.FunctionSpace,
548 ):
549 r"""
550 Tie the dofs on mesh entities rigidly to a point, as the RBE2 element of other codes
551 (a rigid "spider").
553 Each dof of this constraint's space on the entities is a "foot" of a spider whose "body"
554 is a point of the point mesh of `W`. Every component of a foot follows the motion of its
555 body,
557 .. math::
559 u(x) = t + \theta \times (x - x_c),
561 where :math:`x_c` is the coordinate of the dofs of `W` at the point, :math:`t` its
562 translation and :math:`\theta` its rotation, the dofs of `W` at the point. Without
563 rotations, :math:`u(x) = t`. All rotation terms are kept, also where their coefficient is
564 zero, so :meth:`update_rbe2` can follow the motion of the meshes.
566 Args:
567 dim: Topological dimension of the entities
568 entities: Entities (local to the process) whose dofs are tied to spider 0, or a
569 sequence whose entry `k` holds the entities tied to the spider with input index
570 `k` (see :func:`dolfinx_mpc.create_spider_mesh`). An entry may be `None`.
571 W: Space on the spider mesh (:func:`dolfinx_mpc.create_spider_mesh`). Its value size
572 is the geometric dimension, for translations only, or 6 in 3D and 3 in 2D, for
573 translations and rotations. `W` must be the space of another constraint
574 finalized together with this one by :func:`finalize_multipointconstraints`.
576 Note:
577 Collective. Must be called by every process, with the same number of entries in
578 `entities`.
579 """
580 self._raise_if_finalized()
581 per_spider = [entities] if isinstance(entities, numpy.ndarray) else list(entities)
582 dofs = []
583 # The dofs of each spider in turn, collectively, so that a dof on the entities of two
584 # spiders is caught as constrained twice
585 for e in per_spider:
586 e = numpy.zeros(0, dtype=numpy.int32) if e is None else numpy.asarray(e, dtype=numpy.int32)
587 dofs.append(_fem.locate_dofs_topological(self.V, dim, e))
588 self._add_rbe2(dofs, W)
590 def add_rbe2_geometrical(
591 self,
592 locators: Union[
593 Callable[[numpy.ndarray], numpy.ndarray], Sequence[Optional[Callable[[numpy.ndarray], numpy.ndarray]]]
594 ],
595 W: _fem.FunctionSpace,
596 ):
597 r"""
598 Tie the dofs located geometrically rigidly to a point, as the RBE2 element of other codes
599 (a rigid "spider"). See :meth:`add_rbe2_topological` for the relation.
601 Args:
602 locators: Marks the dofs tied to spider 0, given their coordinates, shape
603 `(3, num_points)`, or a sequence whose entry `k` marks the dofs tied to the spider
604 with input index `k`. An entry may be `None`.
605 W: Space on the spider mesh, see :meth:`add_rbe2_topological`
607 Note:
608 Collective. Must be called by every process, with the same number of locators.
609 """
610 self._raise_if_finalized()
611 per_spider = [locators] if callable(locators) else list(locators)
612 dofs = [
613 numpy.zeros(0, dtype=numpy.int32)
614 if locator is None
615 else numpy.asarray(_fem.locate_dofs_geometrical(self.V, locator), dtype=numpy.int32)
616 for locator in per_spider
617 ]
618 self._add_rbe2(dofs, W, self.V.tabulate_dof_coordinates())
620 def update_rbe2(self) -> None:
621 """
622 Recompute the coefficients of every RBE2 constraint from the current coordinates.
624 The feet are at the dof coordinates of the constraint's space, the spiders at those of the
625 space on the spider mesh, both read now. Move the meshes, for instance to the deformed
626 configuration in an updated Lagrangian analysis, then call this to tie the feet to the
627 rigid motion about the new positions. Assemble again afterwards.
629 The constraint must be finalized without a `filter`, which could drop a master whose
630 coefficient becomes nonzero.
632 Note:
633 Collective. Must be called by every process.
634 """
635 self._raise_if_not_finalized()
636 if len(self._rbe2) == 0:
637 raise ValueError("The constraint has no RBE2 constraints")
638 for W, V, block in self._rbe2:
639 dolfinx_mpc.cpp.mpc.update_rbe2(self._cpp_object, V._cpp_object, W._cpp_object, block)
641 def _add_rbe3(self, V: _fem.FunctionSpace, dofs: list[npt.NDArray[numpy.int32]], weights, x):
642 """Record `dofs[k]` as feet of spider `k`, with weights from `weights`."""
643 real = V.mesh.geometry.x.dtype
644 for k, d in enumerate(dofs):
645 if weights is None:
646 w = numpy.ones(len(d), dtype=real)
647 elif callable(weights):
648 w = numpy.asarray(weights(x[d].T), dtype=real).reshape(-1)
649 else:
650 w = numpy.full(len(d), weights, dtype=real)
651 self._rbe3.append((V, d, numpy.full(len(d), k, dtype=numpy.int64), w))
653 def add_rbe3_topological(
654 self,
655 V: _fem.FunctionSpace,
656 dim: int,
657 entities: Union[npt.NDArray[numpy.int32], Sequence[Optional[npt.NDArray[numpy.int32]]]],
658 weights: Union[None, float, Callable[[numpy.ndarray], numpy.ndarray]] = None,
659 ):
660 r"""
661 Tie the dofs of spiders to the motion of the dofs of `V` on mesh entities, as the RBE3
662 element of other codes (a flexible "spider").
664 This constraint is on the space of the spider mesh (:func:`dolfinx_mpc.create_spider_mesh`).
665 Each spider moves with the rigid motion that best fits its "feet", in the weighted
666 least-squares sense,
668 .. math::
670 \min_{t, \theta} \sum_i w_i |u_i - t - \theta \times (x_i - x_c)|^2,
672 where :math:`x_c` is the coordinate of the spider, :math:`t` and :math:`\theta` its
673 translation and rotation, and :math:`u_i` the displacement of foot :math:`i` at
674 :math:`x_i`. Without rotations, :math:`t` is the weighted mean of the feet. Unlike RBE2, the
675 feet keep their stiffness: a load on the spider is spread over them without making them
676 rigid.
678 The feet may be in several spaces, given by one call each. The constraint is built when
679 it is finalized, by :func:`finalize_multipointconstraints` together with the constraints of
680 the spaces of the feet.
682 Args:
683 V: The space of the feet, with one component per dimension
684 dim: Topological dimension of the entities
685 entities: Entities (local to the process) whose dofs are feet of spider 0, or a
686 sequence whose entry `k` holds those of the spider with input index `k`. An entry
687 may be `None`.
688 weights: The weight of each foot: `None` for one, a number, or a function of the
689 coordinates, shape `(3, num_points)`, returning one non-negative weight per foot.
690 Evaluated once, here: :meth:`update_rbe3` keeps the weights.
692 Note:
693 Collective. Must be called by every process, with the same number of entries in
694 `entities`.
695 """
696 self._raise_if_finalized()
697 per_spider = [entities] if isinstance(entities, numpy.ndarray) else list(entities)
698 dofs = []
699 for e in per_spider:
700 e = numpy.zeros(0, dtype=numpy.int32) if e is None else numpy.asarray(e, dtype=numpy.int32)
701 dofs.append(_fem.locate_dofs_topological(V, dim, e))
702 self._add_rbe3(V, dofs, weights, V.tabulate_dof_coordinates() if callable(weights) else None)
704 def add_rbe3_geometrical(
705 self,
706 V: _fem.FunctionSpace,
707 locators: Union[
708 Callable[[numpy.ndarray], numpy.ndarray], Sequence[Optional[Callable[[numpy.ndarray], numpy.ndarray]]]
709 ],
710 weights: Union[None, float, Callable[[numpy.ndarray], numpy.ndarray]] = None,
711 ):
712 """
713 Tie the dofs of spiders to the motion of the dofs of `V` located geometrically, as the
714 RBE3 element of other codes. See :meth:`add_rbe3_topological` for the relation.
716 Args:
717 V: The space of the feet, with one component per dimension
718 locators: Marks the feet of spider 0, given their coordinates, shape
719 `(3, num_points)`, or a sequence whose entry `k` marks the feet of the spider with
720 input index `k`. An entry may be `None`.
721 weights: See :meth:`add_rbe3_topological`
723 Note:
724 Collective. Must be called by every process, with the same number of locators.
725 """
726 self._raise_if_finalized()
727 per_spider = [locators] if callable(locators) else list(locators)
728 dofs = [
729 numpy.zeros(0, dtype=numpy.int32)
730 if locator is None
731 else numpy.asarray(_fem.locate_dofs_geometrical(V, locator), dtype=numpy.int32)
732 for locator in per_spider
733 ]
734 self._add_rbe3(V, dofs, weights, V.tabulate_dof_coordinates() if callable(weights) else None)
736 def _build_rbe3(self, mpcs: Sequence[MultiPointConstraint]) -> None:
737 """Build the RBE3 constraint from the feet recorded, with the blocks of their spaces."""
738 spaces: list[_fem.FunctionSpace] = []
739 for V, *_ in self._rbe3:
740 if not any(V is other for other in spaces):
741 spaces.append(V)
742 blocks = []
743 for V in spaces:
744 matches = [j for j, other in enumerate(mpcs) if other.V is V]
745 if len(matches) != 1:
746 raise ValueError(
747 "The feet of an RBE3 constraint must be in the space of exactly one of the "
748 "constraints finalized together with it"
749 )
750 blocks.append(matches[0])
752 def gather(i):
753 return [numpy.concatenate([r[i] for r in self._rbe3 if r[0] is V]) for V in spaces]
755 data = (self.V, spaces, gather(1), gather(2), gather(3))
756 mpc_data, space = create_rbe3(*data, dtype=self._dtype)
757 self.add_constraint(
758 self.V,
759 mpc_data.slaves,
760 mpc_data.masters,
761 mpc_data.coeffs,
762 mpc_data.owners,
763 mpc_data.offsets,
764 master_blocks=numpy.asarray(blocks, dtype=numpy.int32)[space],
765 )
766 self._rbe3_data = (*data, blocks)
768 def update_rbe3(self) -> None:
769 """
770 Recompute the coefficients of the RBE3 constraint from the current coordinates.
772 The feet are at the dof coordinates of their spaces, the spiders at those of the space of
773 this constraint, both read now. Move the meshes, then call this. Assemble again
774 afterwards.
776 Note:
777 Collective. Must be called by every process.
778 """
779 self._raise_if_not_finalized()
780 if self._rbe3_data is None:
781 raise ValueError("The constraint has no RBE3 constraint")
782 W, spaces, dofs, spiders, weights, blocks = self._rbe3_data
783 dolfinx_mpc.cpp.mpc.update_rbe3(
784 self._cpp_object,
785 W._cpp_object,
786 [V._cpp_object for V in spaces],
787 blocks,
788 [numpy.ascontiguousarray(d, dtype=numpy.int32) for d in dofs],
789 [numpy.ascontiguousarray(k, dtype=numpy.int64) for k in spiders],
790 [numpy.ascontiguousarray(w) for w in weights],
791 )
793 def create_slip_constraint(
794 self,
795 space: _fem.FunctionSpace,
796 facet_marker: Tuple[_mesh.MeshTags, int],
797 v: _fem.Function,
798 bcs: List[_fem.DirichletBC] = [],
799 ):
800 """
801 Create a slip constraint :math:`u \\cdot v=0` over the entities defined in `facet_marker` with the given index.
803 Args:
804 space: Function space (possible sub space) for the current constraint
805 facet_marker: Tuple containomg the mesh tag and marker used to locate degrees of freedom
806 v: Function containing the directional vector to dot your slip condition (most commonly a normal vector)
807 bcs: List of Dirichlet BCs (slip conditions will be ignored on these dofs)
809 Examples:
810 Create constaint :math:`u\\cdot n=0` of all indices in `mt` marked with `i`
812 .. highlight:: python
813 .. code-block:: python
815 V = dolfinx.fem.functionspace(mesh, ("CG", 1))
816 mpc = MultiPointConstraint(V)
817 n = dolfinx.fem.Function(V)
818 mpc.create_slip_constaint(V, (mt, i), n)
820 Create slip constaint for a mixed function space:
822 .. highlight:: python
823 .. code-block:: python
825 cellname = mesh.basix_cell()
826 Ve = basix.ufl.element(basix.ElementFamily.P, cellname , 2, shape=(mesh.geometry.dim,))
827 Qe = basix.ufl.element(basix.ElementFamily.P, cellname , 1)
828 me = basix.ufl.mixed_element([Ve, Qe])
829 W = dolfinx.fem.functionspace(mesh, me)
830 mpc = MultiPointConstraint(W)
831 n_space, _ = W.sub(0).collapse()
832 normal = dolfinx.fem.Function(n_space)
833 mpc.create_slip_constraint(W.sub(0), (mt, i), normal, bcs=[])
835 A slip condition cannot be applied on the same degrees of freedom as a Dirichlet BC, and therefore
836 any Dirichlet bc for the space of the multi point constraint should be supplied.
838 .. highlight:: python
839 .. code-block:: python
841 cellname = mesh.basix_cell()
842 Ve = basix.ufl.element(basix.ElementFamily.P, cellname , 2, shape=(mesh.geometry.dim,))
843 Qe = basix.ufl.element(basix.ElementFamily.P, cellname , 1)
844 me = basix.ufl.mixed_element([Ve, Qe])
845 W = dolfinx.fem.functionspace(mesh, me)
846 mpc = MultiPointConstraint(W)
847 n_space, _ = W.sub(0).collapse()
848 normal = Function(n_space)
849 bc = dolfinx.fem.dirichletbc(inlet_velocity, dofs, W.sub(0))
850 mpc.create_slip_constraint(W.sub(0), (mt, i), normal, bcs=[bc])
851 """
852 bcs = [] if bcs is None else [bc._cpp_object for bc in bcs]
853 if space is self.V:
854 sub_space = False
855 elif self.V.contains(space):
856 sub_space = True
857 else:
858 raise ValueError("Input space has to be a sub space of the MPC space")
859 mpc_data = _cpp_function("create_slip_condition", self._dtype)(
860 space._cpp_object,
861 facet_marker[0]._cpp_object,
862 facet_marker[1],
863 v._cpp_object,
864 bcs,
865 sub_space,
866 )
867 self.add_constraint_from_mpc_data(self.V, mpc_data=mpc_data)
869 def create_general_constraint(
870 self,
871 slave_master_dict: Dict[bytes, Dict[bytes, float]],
872 subspace_slave: Optional[int] = None,
873 subspace_master: Optional[int] = None,
874 *,
875 distance_tol: Optional[float] = None,
876 ):
877 """
878 Args:
879 V: The function space
880 slave_master_dict: Nested dictionary, where the first key is the bit representing the slave dof's
881 coordinate in the mesh. The item of this key is a dictionary, where each key of this dictionary
882 is the bit representation of the master dof's coordinate, and the item the coefficient for
883 the MPC equation.
884 subspace_slave: If using mixed or vector space, and only want to use dofs from a sub space
885 as slave add index here
886 subspace_master: Subspace index for mixed or vector spaces
887 distance_tol: The largest distance between the coordinate of a key and that of its dof.
888 Defaults to `500` machine epsilon of the coordinate type of the mesh.
890 Example:
891 If the dof `D` located at `[d0, d1]` should be constrained to the dofs
892 `E` and `F` at `[e0, e1]` and `[f0, f1]` as :math:`D = \\alpha E + \\beta F`
893 the dictionary should be:
895 .. highlight:: python
896 .. code-block:: python
898 {numpy.array([d0, d1], dtype=mesh.geometry.x.dtype).tobytes():
899 {numpy.array([e0, e1], dtype=mesh.geometry.x.dtype).tobytes(): alpha,
900 numpy.array([f0, f1], dtype=mesh.geometry.x.dtype).tobytes(): beta}}
901 """
902 slaves, masters, coeffs, owners, offsets = create_dictionary_constraint(
903 self.V,
904 slave_master_dict,
905 subspace_slave,
906 subspace_master,
907 dtype=self._dtype,
908 distance_tol=_tolerance(distance_tol, self.V.mesh.geometry.x.dtype),
909 )
910 self.add_constraint(self.V, slaves, masters, coeffs, owners, offsets)
912 def _tolerances(
913 self,
914 distance_tol: Optional[float],
915 coefficient_tol: Optional[float],
916 tol: Union[_float_classes, float, None, _Unset] = _UNSET,
917 eps2: Optional[float] = None,
918 ) -> tuple[float, float]:
919 """The distance and coefficient tolerance, by default `500` machine epsilon of the coordinate
920 type of the mesh and of the real type of the constraint.
922 The deprecated `tol` of the periodic constraints sets both, or with `None` keeps every
923 master. The deprecated `eps2` of the contact constraints is a squared distance.
924 """
925 if not isinstance(tol, _Unset):
926 _deprecated("tol", "`distance_tol` and `coefficient_tol`", stacklevel=4)
927 if tol is None:
928 coefficient_tol = 0.0 if coefficient_tol is None else coefficient_tol
929 else:
930 distance_tol = tol if distance_tol is None else distance_tol
931 coefficient_tol = tol if coefficient_tol is None else coefficient_tol
932 if eps2 is not None:
933 _deprecated("eps2", "`distance_tol`, a distance rather than a squared distance,", stacklevel=4)
934 distance_tol = float(numpy.sqrt(eps2)) if distance_tol is None else distance_tol
935 return (
936 _tolerance(distance_tol, self.V.mesh.geometry.x.dtype),
937 _tolerance(coefficient_tol, self._dtype),
938 )
940 def create_contact_slip_condition(
941 self,
942 meshtags: _mesh.MeshTags,
943 slave_marker: int,
944 master_marker: int,
945 normal: _fem.Function,
946 eps2: Optional[float] = None,
947 num_threads: Optional[int] = 1,
948 *,
949 distance_tol: Optional[float] = None,
950 coefficient_tol: Optional[float] = None,
951 ):
952 """
953 Create a slip condition between two sets of facets marker with individual markers.
954 The interfaces should be within machine precision of eachother, but the vertices does not need to align.
955 The condition created is :math:`u_s \\cdot normal_s = u_m \\cdot normal_m` where `s` is the
956 restriction to the slave facets, `m` to the master facets.
958 Args:
959 meshtags: The meshtags of the set of facets to tie together
960 slave_marker: The marker of the slave facets
961 master_marker: The marker of the master facets
962 normal: The function used in the dot-product of the constraint
963 eps2: Deprecated, use `distance_tol`, which is `sqrt(eps2)`.
964 num_threads: The number of threads to use for certain operations
965 distance_tol: The largest distance from a slave point to a master cell for the point to
966 be in the cell, and the padding of the bounding boxes of the cells. Defaults to `500`
967 machine epsilon of the coordinate type of the mesh.
968 coefficient_tol: A master whose coefficient is below `coefficient_tol` times the largest of
969 its slave is dropped. `0` keeps every master. Defaults to `500`
970 machine epsilon of the real type of the constraint.
971 """
972 mpc_data = _cpp_function("create_contact_slip_condition", self._dtype)(
973 self.V._cpp_object,
974 meshtags._cpp_object,
975 slave_marker,
976 master_marker,
977 normal._cpp_object,
978 *self._tolerances(distance_tol, coefficient_tol, eps2=eps2),
979 num_threads,
980 )
981 self.add_constraint_from_mpc_data(self.V, mpc_data)
983 def create_contact_inelastic_condition(
984 self,
985 meshtags: _cpp.mesh.MeshTags_int32,
986 slave_marker: int,
987 master_marker: int,
988 eps2: Optional[float] = None,
989 allow_missing_masters: bool = False,
990 num_threads: Optional[int] = 1,
991 *,
992 distance_tol: Optional[float] = None,
993 coefficient_tol: Optional[float] = None,
994 ):
995 """
996 Create a contact inelastic condition between two sets of facets marker with individual markers.
997 The interfaces should be within machine precision of eachother, but the vertices does not need to align.
998 The condition created is :math:`u_s = u_m` where `s` is the restriction to the
999 slave facets, `m` to the master facets.
1001 Args:
1002 meshtags: The meshtags of the set of facets to tie together
1003 slave_marker: The marker of the slave facets
1004 master_marker: The marker of the master facets
1005 eps2: Deprecated, use `distance_tol`, which is `sqrt(eps2)`.
1006 allow_missing_masters: If true, the function will not throw an error if a degree of freedom
1007 in the closure of the master entities does not have a corresponding set of slave degree
1008 of freedom.
1009 num_threads: The number of threads to use for certain operations
1010 distance_tol: The largest distance from a slave point to a master cell for the point to
1011 be in the cell, and the padding of the bounding boxes of the cells. Defaults to `500`
1012 machine epsilon of the coordinate type of the mesh.
1013 coefficient_tol: A master whose coefficient is below `coefficient_tol` times the largest of
1014 its slave is dropped. `0` keeps every master. Defaults to `500`
1015 machine epsilon of the real type of the constraint.
1016 """
1017 mpc_data = _cpp_function("create_contact_inelastic_condition", self._dtype)(
1018 self.V._cpp_object,
1019 meshtags._cpp_object,
1020 slave_marker,
1021 master_marker,
1022 *self._tolerances(distance_tol, coefficient_tol, eps2=eps2),
1023 allow_missing_masters,
1024 num_threads,
1025 )
1026 self.add_constraint_from_mpc_data(self.V, mpc_data)
1028 @property
1029 def is_slave(self) -> numpy.ndarray:
1030 """
1031 Returns a vector of integers where the ith entry indicates if a degree of freedom (local to process) is a slave.
1032 """
1033 self._raise_if_not_finalized()
1034 return self._cpp_object.is_slave
1036 @property
1037 def slaves(self):
1038 """
1039 Returns the degrees of freedom for all slaves local to process
1040 """
1041 self._raise_if_not_finalized()
1042 return self._cpp_object.slaves
1044 @property
1045 def masters(self) -> _cpp.graph.AdjacencyList_int32:
1046 """
1047 Returns an adjacency-list whose ith node corresponds to
1048 a degree of freedom (local to process), and links the corresponding master dofs (local to process).
1050 Examples:
1052 .. highlight:: python
1053 .. code-block:: python
1055 masters = mpc.masters
1056 masters_of_dof_i = masters.links(i)
1057 """
1058 self._raise_if_not_finalized()
1059 return self._cpp_object.masters
1061 def coefficients(self) -> _float_array_types:
1062 """
1063 Returns a vector containing the coefficients for the constraint, and the corresponding offsets
1064 for the ith degree of freedom.
1066 Examples:
1068 .. highlight:: python
1069 .. code-block:: python
1071 coeffs, offsets = mpc.coefficients()
1072 coeffs_of_slave_i = coeffs[offsets[i]:offsets[i+1]]
1073 """
1074 self._raise_if_not_finalized()
1075 return self._cpp_object.coefficients()
1077 def all_coefficients(self) -> Tuple[_float_array_types, npt.NDArray[numpy.int32]]:
1078 """
1079 Returns the coefficients of all masters, including those eliminated by a Dirichlet condition,
1080 in the order supplied before :func:`finalize`, and the offsets for the ith degree of freedom.
1081 This is the layout taken by :func:`update_coefficients`. The corresponding masters are given
1082 by :func:`all_masters`.
1084 Examples:
1086 .. highlight:: python
1087 .. code-block:: python
1089 coeffs, offsets = mpc.all_coefficients()
1090 coeffs_of_slave_i = coeffs[offsets[i]:offsets[i+1]]
1091 """
1092 self._raise_if_not_finalized()
1093 return self._cpp_object.all_coefficients()
1095 def all_masters(self) -> npt.NDArray[numpy.int32]:
1096 """
1097 Returns the masters (local index in :attr:`function_space`) in the layout of
1098 :func:`all_coefficients`.
1099 """
1100 self._raise_if_not_finalized()
1101 return self._cpp_object.all_masters()
1103 def update_coefficients(self, coeffs: _float_array_types) -> None:
1104 """
1105 Replace the coefficient of every master, including masters eliminated by a Dirichlet
1106 condition, and recompute the constraint offset :math:`g`.
1108 The masters are fixed at creation. A master dropped by `coefficient_tol` or by the `filter`
1109 of :func:`finalize` cannot be given a coefficient, so create the constraint with
1110 `coefficient_tol=0` and no filter if the coefficients are to be changed.
1112 Args:
1113 coeffs: The new coefficients, in the layout of :func:`all_coefficients`, for all degrees
1114 of freedom local to the process (owned and ghost).
1116 Note:
1117 Collective. Must be called by every process.
1118 """
1119 self._raise_if_not_finalized()
1120 self._cpp_object.update_coefficients(numpy.ascontiguousarray(coeffs, dtype=self._dtype))
1122 def scale_coefficients(
1123 self,
1124 scale: Union[_float_classes, float, complex, ufl.core.expr.Expr, _fem.Expression],
1125 ) -> None:
1126 """
1127 Multiply the coefficients of all masters of each slave :math:`s` by a factor
1128 :math:`f_s`, and recompute the constraint offset :math:`g`. For a periodic constraint
1129 :math:`u(x_s) = f_s u(relation(x_s))`, which for instance gives a Floquet-Bloch condition
1130 with :math:`f=e^{i k\\cdot L}`.
1132 The factors are stored in a function in the space of the constraint, and :math:`f_s` is
1133 the degree of freedom :math:`s` of that function: the value at the slave for a Lagrange
1134 space, the corresponding moment for e.g. a Nédélec space.
1136 Repeated calls compound. Masters eliminated by a Dirichlet condition are scaled as well,
1137 the user supplied `rhs_coeffs` are not.
1139 Args:
1140 scale: A scalar, a :class:`dolfinx.fem.Function` in the constraint's space (copied
1141 by interpolation), a UFL expression, compiled into a :class:`dolfinx.fem.Expression`
1142 at the interpolation points of the space, or such a compiled expression. Pass a
1143 compiled expression to avoid recompilation when the factor is updated through
1144 :class:`dolfinx.fem.Constant`'s in it.
1146 Note:
1147 Collective. Must be called by every process.
1148 """
1149 self._raise_if_not_finalized()
1150 if self._scale_function is None:
1151 self._scale_function = _fem.Function(self.V, dtype=self._dtype)
1152 f = self._scale_function
1153 if isinstance(scale, (_fem.Expression, _fem.Function)):
1154 f.interpolate(scale)
1155 elif isinstance(scale, ufl.core.expr.Expr):
1156 f.interpolate(_fem.Expression(scale, self.V.element.interpolation_points, dtype=self._dtype))
1157 else:
1158 f.x.array[:] = scale
1159 f.x.scatter_forward()
1160 # The extended index map appends master ghosts after the ghosts of the input space
1161 num_dofs_local = len(self._cpp_object.is_slave)
1162 self._cpp_object.scale_coefficients(f.x.array[:num_dofs_local])
1164 @property
1165 def num_local_slaves(self):
1166 """
1167 Return the number of slaves owned by the current process.
1168 """
1169 self._raise_if_not_finalized()
1170 return self._cpp_object.num_local_slaves
1172 @property
1173 def cell_to_slaves(self):
1174 """
1175 Returns an `dolfinx.cpp.graph.AdjacencyList_int32` whose ith node corresponds to
1176 the ith cell (local to process), and links the corresponding slave degrees of
1177 freedom in the cell (local to process).
1179 Examples:
1181 .. highlight:: python
1182 .. code-block:: python
1184 cell_to_slaves = mpc.cell_to_slaves()
1185 slaves_in_cell_i = cell_to_slaves.links(i)
1186 """
1187 self._raise_if_not_finalized()
1188 return self._cpp_object.cell_to_slaves
1190 @property
1191 def function_space(self):
1192 """
1193 Return the function space for the multi-point constraint with the updated index map
1194 """
1195 self._raise_if_not_finalized()
1196 return self.V
1198 @property
1199 def input_space(self) -> _fem.FunctionSpace:
1200 """
1201 The function space the constraint was created with.
1203 Forms and Dirichlet conditions are stated on this space, while functions holding a
1204 solution live in :attr:`function_space`, its extension by the masters of the constraint.
1205 For a system of several blocks, ``[mpc.input_space for mpc in mpcs]`` gives the spaces in
1206 the order of the blocks, for instance for a :class:`ufl.MixedFunctionSpace`.
1207 """
1208 return self._input_space
1210 def backsubstitution(self, u: Union[_fem.Function, Sequence[_fem.Function], _PETSc.Vec]) -> None: # type: ignore
1211 """
1212 For a Function, impose the multi-point constraint by backsubstiution.
1213 This function is used after solving the reduced problem to obtain the values
1214 at the slave degrees of freedom
1216 .. note::
1217 It is the users responsibility to destroy the PETSc vector
1219 Args:
1220 u: The input function. For a constraint with masters in another block, the function
1221 of every block, in the order given to :func:`finalize_multipointconstraints`;
1222 only the function of this constraint's block is changed. The ghosts of the
1223 functions holding masters must be up to date.
1224 """
1225 self._raise_if_not_finalized()
1226 if isinstance(u, Sequence):
1227 self._cpp_object.backsubstitution([u_k.x.array for u_k in u]) # type: ignore
1228 u[self._cpp_object.block].x.scatter_forward()
1229 return
1230 try:
1231 self._cpp_object.backsubstitution(u.x.array) # type: ignore
1232 assert isinstance(u, _fem.Function)
1233 u.x.scatter_forward()
1234 except AttributeError:
1235 assert isinstance(u, _PETSc.Vec)
1236 with u.localForm() as vector_local:
1237 self._cpp_object.backsubstitution(vector_local.array_w)
1238 u.ghostUpdate(addv=_PETSc.InsertMode.INSERT, mode=_PETSc.ScatterMode.FORWARD) # type: ignore
1240 def homogenize(self, u: _fem.Function) -> None:
1241 """
1242 For a vector, homogenize (set to zero) the vector components at the multi-point
1243 constraint slave DoF indices. This is particularly useful for nonlinear problems.
1245 Args:
1246 u: The input vector
1247 """
1248 self._cpp_object.homogenize(u.x.array)
1249 u.x.scatter_forward()
1251 def _raise_if_finalized(self):
1252 """
1253 Raise if the multi point constraint has already been finalized
1254 """
1255 if self.finalized:
1256 raise RuntimeError("MultiPointConstraint has already been finalized")
1258 def _raise_if_not_finalized(self):
1259 """
1260 Raise if the multi point constraint has not yet been finalized
1261 """
1262 if not self.finalized:
1263 raise RuntimeError("MultiPointConstraint has not been finalized")
1266def finalize_multipointconstraints(
1267 mpcs: Sequence[MultiPointConstraint], filter: Optional[numpy.floating] = None
1268) -> None:
1269 """
1270 Finalize the multi point constraints of several function spaces together.
1272 Entry ``k`` of ``mpcs`` constrains its own function space, for instance the ``k``-th block of a
1273 :class:`ufl.MixedFunctionSpace`. Each is finalized as by :meth:`MultiPointConstraint.finalize`,
1274 but the checks that need communication are reduced once for all of them, and the meshes of the
1275 spaces may be distinct, as long as they live on congruent communicators.
1277 Args:
1278 mpcs: The constraints to finalize. None may be finalized already, and they must all use the
1279 same ``dtype``.
1280 filter: See :meth:`MultiPointConstraint.finalize`. Applied to every constraint.
1282 Raises:
1283 ValueError: If the input is inconsistent, or if a dof is both a slave and constrained by a
1284 Dirichlet condition, a master is also a slave, or the meshes are on communicators
1285 of different size or rank order. Raised on every process.
1287 Note:
1288 Collective. Must be called by every process, with the constraints in the same order.
1289 """
1290 mpcs = list(mpcs)
1291 if len(mpcs) == 0:
1292 raise ValueError("At least one constraint is required")
1293 if len({id(mpc) for mpc in mpcs}) != len(mpcs):
1294 raise ValueError("The same constraint was given more than once")
1295 for mpc in mpcs:
1296 mpc._raise_if_finalized()
1297 dtype = numpy.dtype(mpcs[0]._dtype)
1298 if any(numpy.dtype(mpc._dtype) != dtype for mpc in mpcs):
1299 raise ValueError("All constraints must have the same dtype")
1300 if dtype.type not in (numpy.float32, numpy.float64, numpy.complex64, numpy.complex128):
1301 raise ValueError(f"Unsupported dtype {dtype} for coefficients")
1303 # An RBE3 constraint is built now, when the block of each space of its feet is known
1304 for mpc in mpcs:
1305 if len(mpc._rbe3) > 0:
1306 mpc._build_rbe3(mpcs)
1308 rhs_coeffs = []
1309 for mpc in mpcs:
1310 if mpc._rhs_coeffs is None:
1311 rhs_coeffs.append(numpy.zeros(0, dtype=dtype))
1312 else:
1313 num_dofs_local = mpc.V.dofmap.index_map_bs * (
1314 mpc.V.dofmap.index_map.size_local + mpc.V.dofmap.index_map.num_ghosts
1315 )
1316 rhs_coeffs.append(mpc._rhs_coeffs.x.array[:num_dofs_local].astype(dtype))
1318 # The block of each master: -2 marks the constraint's own block, -1 the block of the space
1319 # recorded with it, anything else a block given directly. Every process records the same chunks
1320 # with the same spaces, so a space that is not one of the blocks raises everywhere.
1321 master_blocks = []
1322 for k, mpc in enumerate(mpcs):
1323 if all(space is None and (blocks == -2).all() for blocks, space in mpc._master_spaces):
1324 master_blocks.append(numpy.zeros(0, dtype=numpy.int32))
1325 continue
1326 resolved = []
1327 for blocks, space in mpc._master_spaces:
1328 blocks = blocks.copy()
1329 blocks[blocks == -2] = k
1330 if space is not None:
1331 matches = [j for j, other in enumerate(mpcs) if other.V is space]
1332 if len(matches) == 0:
1333 matches = [j for j, other in enumerate(mpcs) if other.V.contains(space)]
1334 if len(matches) != 1:
1335 raise ValueError(
1336 "The master space of a constraint must be the function space, or a subspace of "
1337 "the function space, of exactly one of the constraints finalized together with it"
1338 )
1339 blocks[blocks == -1] = matches[0]
1340 resolved.append(blocks)
1341 master_blocks.append(numpy.concatenate(resolved) if resolved else numpy.zeros(0, dtype=numpy.int32))
1343 # Raises ValueError (as the C++ throws std::invalid_argument), identically on every process
1344 cpp_objects = dolfinx_mpc.cpp.mpc.create_multipointconstraints(
1345 [mpc.V._cpp_object for mpc in mpcs],
1346 [mpc._slaves for mpc in mpcs],
1347 [mpc._masters for mpc in mpcs],
1348 [mpc._coeffs.astype(dtype) for mpc in mpcs],
1349 [mpc._owners for mpc in mpcs],
1350 [mpc._offsets for mpc in mpcs],
1351 rhs_coeffs,
1352 [[bc._cpp_object for bc in mpc._bcs] for mpc in mpcs],
1353 master_blocks,
1354 filter,
1355 )
1357 # The block of each space on a spider mesh, matched before the spaces are replaced
1358 for mpc in mpcs:
1359 for entry in mpc._rbe2:
1360 entry[2] = next(j for j, other in enumerate(mpcs) if other.V is entry[0])
1362 for mpc, cpp_object in zip(mpcs, cpp_objects):
1363 mpc._cpp_object = cpp_object
1364 # Replace function space
1365 mpc.V = _fem.FunctionSpace(mpc.V.mesh, mpc.V.ufl_element(), cpp_object.function_space)
1366 mpc.finalized = True
1367 # Delete variables that are no longer required
1368 del (mpc._slaves, mpc._masters, mpc._coeffs, mpc._owners, mpc._offsets, mpc._master_spaces)