Coverage for python/src/dolfinx_mpc/multipointconstraint.py: 82%
175 statements
« prev ^ index » next coverage.py v7.10.1, created at 2025-07-29 22:03 +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: numpy.floating
106 __slots__ = tuple(__annotations__)
108 def __init__(self, V: _fem.FunctionSpace, dtype: numpy.floating = 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.append(self._coeffs, coeffs)
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),
233 tol: _float_classes = 500 * numpy.finfo(default_real_type).eps,
234 ):
235 """
236 Create periodic condition for all closure dofs of on all entities in `meshtag` with value `tag`.
237 :math:`u(x_i) = scale * u(relation(x_i))` for all of :math:`x_i` on marked entities.
239 Args:
240 V: The function space to assign the condition to. Should either be the space of the MPC or a sub space.
241 meshtag: MeshTag for entity to apply the periodic condition on
242 tag: Tag indicating which entities should be slaves
243 relation: Lambda-function describing the geometrical relation
244 bcs: Dirichlet boundary conditions for the problem (Periodic constraints will be ignored for these dofs)
245 scale: Float for scaling bc
246 tol: Tolerance for adding scaled basis values to MPC. Any contribution that is less than this value
247 is ignored. The tolerance is also added as padding for the bounding box trees and corresponding
248 collision searches to determine periodic degrees of freedom.
249 """
250 bcs_ = [bc._cpp_object for bc in bcs]
251 if isinstance(scale, numpy.generic): # nanobind conversion of numpy dtypes to general Python types
252 scale = scale.item() # type: ignore
253 if V is self.V:
254 mpc_data = dolfinx_mpc.cpp.mpc.create_periodic_constraint_topological(
255 self.V._cpp_object, meshtag._cpp_object, tag, relation, bcs_, scale, False, float(tol)
256 )
257 elif self.V.contains(V):
258 mpc_data = dolfinx_mpc.cpp.mpc.create_periodic_constraint_topological(
259 V._cpp_object, meshtag._cpp_object, tag, relation, bcs_, scale, True, float(tol)
260 )
261 else:
262 raise RuntimeError("The input space has to be a sub space (or the full space) of the MPC")
263 self.add_constraint_from_mpc_data(self.V, mpc_data=mpc_data)
265 def create_periodic_constraint_geometrical(
266 self,
267 V: _fem.FunctionSpace,
268 indicator: Callable[[numpy.ndarray], numpy.ndarray],
269 relation: Callable[[numpy.ndarray], numpy.ndarray],
270 bcs: List[_fem.DirichletBC],
271 scale: _float_classes = default_scalar_type(1.0),
272 tol: _float_classes = 500 * numpy.finfo(default_real_type).eps,
273 ):
274 """
275 Create a periodic condition for all degrees of freedom whose physical location satisfies
276 :math:`indicator(x_i)==True`, i.e.
277 :math:`u(x_i) = scale * u(relation(x_i))` for all :math:`x_i`
279 Args:
280 V: The function space to assign the condition to. Should either be the space of the MPC or a sub space.
281 indicator: Lambda-function to locate degrees of freedom that should be slaves
282 relation: Lambda-function describing the geometrical relation to master dofs
283 bcs: Dirichlet boundary conditions for the problem
284 (Periodic constraints will be ignored for these dofs)
285 scale: Float for scaling bc
286 tol: Tolerance for adding scaled basis values to MPC. Any contribution that is less than this value
287 is ignored. The tolerance is also added as padding for the bounding box trees and corresponding
288 collision searches to determine periodic degrees of freedom.
289 """
290 if isinstance(scale, numpy.generic): # nanobind conversion of numpy dtypes to general Python types
291 scale = scale.item() # type: ignore
292 bcs = [] if bcs is None else [bc._cpp_object for bc in bcs]
293 if V is self.V:
294 mpc_data = dolfinx_mpc.cpp.mpc.create_periodic_constraint_geometrical(
295 self.V._cpp_object, indicator, relation, bcs, scale, False, float(tol)
296 )
297 elif self.V.contains(V):
298 mpc_data = dolfinx_mpc.cpp.mpc.create_periodic_constraint_geometrical(
299 V._cpp_object, indicator, relation, bcs, scale, True, float(tol)
300 )
301 else:
302 raise RuntimeError("The input space has to be a sub space (or the full space) of the MPC")
303 self.add_constraint_from_mpc_data(self.V, mpc_data=mpc_data)
305 def create_slip_constraint(
306 self,
307 space: _fem.FunctionSpace,
308 facet_marker: Tuple[_mesh.MeshTags, int],
309 v: _fem.Function,
310 bcs: List[_fem.DirichletBC] = [],
311 ):
312 """
313 Create a slip constraint :math:`u \\cdot v=0` over the entities defined in `facet_marker` with the given index.
315 Args:
316 space: Function space (possible sub space) for the current constraint
317 facet_marker: Tuple containomg the mesh tag and marker used to locate degrees of freedom
318 v: Function containing the directional vector to dot your slip condition (most commonly a normal vector)
319 bcs: List of Dirichlet BCs (slip conditions will be ignored on these dofs)
321 Examples:
322 Create constaint :math:`u\\cdot n=0` of all indices in `mt` marked with `i`
324 .. highlight:: python
325 .. code-block:: python
327 V = dolfinx.fem.functionspace(mesh, ("CG", 1))
328 mpc = MultiPointConstraint(V)
329 n = dolfinx.fem.Function(V)
330 mpc.create_slip_constaint(V, (mt, i), n)
332 Create slip constaint for a mixed function space:
334 .. highlight:: python
335 .. code-block:: python
337 cellname = mesh.ufl_cell().cellname()
338 Ve = basix.ufl.element(basix.ElementFamily.P, cellname , 2, shape=(mesh.geometry.dim,))
339 Qe = basix.ufl.element(basix.ElementFamily.P, cellname , 1)
340 me = basix.ufl.mixed_element([Ve, Qe])
341 W = dolfinx.fem.functionspace(mesh, me)
342 mpc = MultiPointConstraint(W)
343 n_space, _ = W.sub(0).collapse()
344 normal = dolfinx.fem.Function(n_space)
345 mpc.create_slip_constraint(W.sub(0), (mt, i), normal, bcs=[])
347 A slip condition cannot be applied on the same degrees of freedom as a Dirichlet BC, and therefore
348 any Dirichlet bc for the space of the multi point constraint should be supplied.
350 .. highlight:: python
351 .. code-block:: python
353 cellname = mesh.ufl_cell().cellname()
354 Ve = basix.ufl.element(basix.ElementFamily.P, cellname , 2, shape=(mesh.geometry.dim,))
355 Qe = basix.ufl.element(basix.ElementFamily.P, cellname , 1)
356 me = basix.ufl.mixed_element([Ve, Qe])
357 W = dolfinx.fem.functionspace(mesh, me)
358 mpc = MultiPointConstraint(W)
359 n_space, _ = W.sub(0).collapse()
360 normal = Function(n_space)
361 bc = dolfinx.fem.dirichletbc(inlet_velocity, dofs, W.sub(0))
362 mpc.create_slip_constraint(W.sub(0), (mt, i), normal, bcs=[bc])
363 """
364 bcs = [] if bcs is None else [bc._cpp_object for bc in bcs]
365 if space is self.V:
366 sub_space = False
367 elif self.V.contains(space):
368 sub_space = True
369 else:
370 raise ValueError("Input space has to be a sub space of the MPC space")
371 mpc_data = dolfinx_mpc.cpp.mpc.create_slip_condition(
372 space._cpp_object,
373 facet_marker[0]._cpp_object,
374 facet_marker[1],
375 v._cpp_object,
376 bcs,
377 sub_space,
378 )
379 self.add_constraint_from_mpc_data(self.V, mpc_data=mpc_data)
381 def create_general_constraint(
382 self,
383 slave_master_dict: Dict[bytes, Dict[bytes, float]],
384 subspace_slave: Optional[int] = None,
385 subspace_master: Optional[int] = None,
386 ):
387 """
388 Args:
389 V: The function space
390 slave_master_dict: Nested dictionary, where the first key is the bit representing the slave dof's
391 coordinate in the mesh. The item of this key is a dictionary, where each key of this dictionary
392 is the bit representation of the master dof's coordinate, and the item the coefficient for
393 the MPC equation.
394 subspace_slave: If using mixed or vector space, and only want to use dofs from a sub space
395 as slave add index here
396 subspace_master: Subspace index for mixed or vector spaces
398 Example:
399 If the dof `D` located at `[d0, d1]` should be constrained to the dofs
400 `E` and `F` at `[e0, e1]` and `[f0, f1]` as :math:`D = \\alpha E + \\beta F`
401 the dictionary should be:
403 .. highlight:: python
404 .. code-block:: python
406 {numpy.array([d0, d1], dtype=mesh.geometry.x.dtype).tobytes():
407 {numpy.array([e0, e1], dtype=mesh.geometry.x.dtype).tobytes(): alpha,
408 numpy.array([f0, f1], dtype=mesh.geometry.x.dtype).tobytes(): beta}}
409 """
410 slaves, masters, coeffs, owners, offsets = create_dictionary_constraint(
411 self.V, slave_master_dict, subspace_slave, subspace_master
412 )
413 self.add_constraint(self.V, slaves, masters, coeffs, owners, offsets)
415 def create_contact_slip_condition(
416 self,
417 meshtags: _mesh.MeshTags,
418 slave_marker: int,
419 master_marker: int,
420 normal: _fem.Function,
421 eps2: float = 1e-20,
422 ):
423 """
424 Create a slip condition between two sets of facets marker with individual markers.
425 The interfaces should be within machine precision of eachother, but the vertices does not need to align.
426 The condition created is :math:`u_s \\cdot normal_s = u_m \\cdot normal_m` where `s` is the
427 restriction to the slave facets, `m` to the master facets.
429 Args:
430 meshtags: The meshtags of the set of facets to tie together
431 slave_marker: The marker of the slave facets
432 master_marker: The marker of the master facets
433 normal: The function used in the dot-product of the constraint
434 eps2: The tolerance for the squared distance between cells to be considered as a collision
435 """
436 if isinstance(eps2, numpy.generic): # nanobind conversion of numpy dtypes to general Python types
437 eps2 = eps2.item() # type: ignore
438 mpc_data = dolfinx_mpc.cpp.mpc.create_contact_slip_condition(
439 self.V._cpp_object,
440 meshtags._cpp_object,
441 slave_marker,
442 master_marker,
443 normal._cpp_object,
444 eps2,
445 )
446 self.add_constraint_from_mpc_data(self.V, mpc_data)
448 def create_contact_inelastic_condition(
449 self,
450 meshtags: _cpp.mesh.MeshTags_int32,
451 slave_marker: int,
452 master_marker: int,
453 eps2: float = 1e-20,
454 allow_missing_masters: bool = False,
455 ):
456 """
457 Create a contact inelastic condition between two sets of facets marker with individual markers.
458 The interfaces should be within machine precision of eachother, but the vertices does not need to align.
459 The condition created is :math:`u_s = u_m` where `s` is the restriction to the
460 slave facets, `m` to the master facets.
462 Args:
463 meshtags: The meshtags of the set of facets to tie together
464 slave_marker: The marker of the slave facets
465 master_marker: The marker of the master facets
466 eps2: The tolerance for the squared distance between cells to be considered as a collision
467 allow_missing_masters: If true, the function will not throw an error if a degree of freedom
468 in the closure of the master entities does not have a corresponding set of slave degree
469 of freedom.
470 """
471 if isinstance(eps2, numpy.generic): # nanobind conversion of numpy dtypes to general Python types
472 eps2 = eps2.item() # type: ignore
473 mpc_data = dolfinx_mpc.cpp.mpc.create_contact_inelastic_condition(
474 self.V._cpp_object, meshtags._cpp_object, slave_marker, master_marker, eps2, allow_missing_masters
475 )
476 self.add_constraint_from_mpc_data(self.V, mpc_data)
478 @property
479 def is_slave(self) -> numpy.ndarray:
480 """
481 Returns a vector of integers where the ith entry indicates if a degree of freedom (local to process) is a slave.
482 """
483 self._not_finalized()
484 return self._cpp_object.is_slave
486 @property
487 def slaves(self):
488 """
489 Returns the degrees of freedom for all slaves local to process
490 """
491 self._not_finalized()
492 return self._cpp_object.slaves
494 @property
495 def masters(self) -> _cpp.graph.AdjacencyList_int32:
496 """
497 Returns an adjacency-list whose ith node corresponds to
498 a degree of freedom (local to process), and links the corresponding master dofs (local to process).
500 Examples:
502 .. highlight:: python
503 .. code-block:: python
505 masters = mpc.masters
506 masters_of_dof_i = masters.links(i)
507 """
508 self._not_finalized()
509 return self._cpp_object.masters
511 def coefficients(self) -> _float_array_types:
512 """
513 Returns a vector containing the coefficients for the constraint, and the corresponding offsets
514 for the ith degree of freedom.
516 Examples:
518 .. highlight:: python
519 .. code-block:: python
521 coeffs, offsets = mpc.coefficients()
522 coeffs_of_slave_i = coeffs[offsets[i]:offsets[i+1]]
523 """
524 self._not_finalized()
525 return self._cpp_object.coefficients()
527 @property
528 def num_local_slaves(self):
529 """
530 Return the number of slaves owned by the current process.
531 """
532 self._not_finalized()
533 return self._cpp_object.num_local_slaves
535 @property
536 def cell_to_slaves(self):
537 """
538 Returns an `dolfinx.cpp.graph.AdjacencyList_int32` whose ith node corresponds to
539 the ith cell (local to process), and links the corresponding slave degrees of
540 freedom in the cell (local to process).
542 Examples:
544 .. highlight:: python
545 .. code-block:: python
547 cell_to_slaves = mpc.cell_to_slaves()
548 slaves_in_cell_i = cell_to_slaves.links(i)
549 """
550 self._not_finalized()
551 return self._cpp_object.cell_to_slaves
553 @property
554 def function_space(self):
555 """
556 Return the function space for the multi-point constraint with the updated index map
557 """
558 self._not_finalized()
559 return self.V
561 def backsubstitution(self, u: Union[_fem.Function, _PETSc.Vec]) -> None: # type: ignore
562 """
563 For a Function, impose the multi-point constraint by backsubstiution.
564 This function is used after solving the reduced problem to obtain the values
565 at the slave degrees of freedom
567 .. note::
568 It is the users responsibility to destroy the PETSc vector
570 Args:
571 u: The input function
572 """
573 try:
574 self._cpp_object.backsubstitution(u.x.array)
575 u.x.scatter_forward()
576 except AttributeError:
577 with u.localForm() as vector_local:
578 self._cpp_object.backsubstitution(vector_local.array_w)
579 u.ghostUpdate(addv=_PETSc.InsertMode.INSERT, mode=_PETSc.ScatterMode.FORWARD) # type: ignore
581 def homogenize(self, u: _fem.Function) -> None:
582 """
583 For a vector, homogenize (set to zero) the vector components at the multi-point
584 constraint slave DoF indices. This is particularly useful for nonlinear problems.
586 Args:
587 u: The input vector
588 """
589 self._cpp_object.homogenize(u.x.array)
590 u.x.scatter_forward()
592 def _already_finalized(self):
593 """
594 Check if we have already finalized the multi point constraint
595 """
596 if self.finalized:
597 raise RuntimeError("MultiPointConstraint has already been finalized")
599 def _not_finalized(self):
600 """
601 Check if we have finalized the multi point constraint
602 """
603 if not self.finalized:
604 raise RuntimeError("MultiPointConstraint has not been finalized")