Boundary conditions and constraints
Basic boundary conditions
There is 3 basic types of boundary conditions in fedoo:
Dirichlet boundary conditions: Value specified for a degree of freedom (for instance displacement).
Neumann boundary conditions: Value specified for the dual value of a degree of freedom (for instance force)
Multi-point constraints (MPC): linear coupling between degrees of freedom. MPC are not properly speaking boundary conditions. In fedoo they are considered as part of boundary conditions, because they often used to specify complexe boundary conditions, but they are more general that simple boundary conditions.
Each time a problem is created, a list of boundary conditions is associated to the problem in the “bc” attribute of the problem instance.
problem.bc is then a fedoo.ListBC Object.
|
List of boundary conditions. |
Dirichlet and Neumann Boundary conditions
The recommanded way to apply some boundary conditions is to use the method
fedoo.ListBC.add(), directly from problem attribute
fedoo.Problem.bc.
This method build a fedoo.core.BoundaryCondition object and add it to
the problem boundary conditions.
|
Basic boundary conditions. |
Here are some examples to add some Dirichlet or Neumann Boundary conditions for a given problem pb:
Displacement along the X axis blocked (=0) for the nodes 0, 1 and 8
>>> pb.bc.add('Dirichlet',[0,1,8], 'DispX', 0)
Displacement along the all axis blocked (=0) for the nodes defined by the node sets ‘node_sets’ (the given node nets is supposed to be present in the associated mesh. If note, it will raise an error).
>>> pb.bc.add('Dirichlet','node_sets', 'Disp', 0)
External force on the node 112 with Fx = 5, Fy = 10 and Fz = 0.
>>> pb.bc.add('Neumann', [112], 'Disp', [5, 10, 0])
Enforce a nul displacement on nodes 1, 2 and 3, and a rotation around X = 0 for node 1, 0.1 for node 2 and 0.2 for node 3.
>>> pb.bc.add('Dirichlet',
>>> [1,2,3],
>>> ['Disp','RotX'],
>>> [0, np.array([0,0.1,0.2])]
>>> )
Instead of adding directly the boundary conditions to the problem, it may be convenient to create boundary conditions and add them afterwards
>>> my_bc = fd.BoundaryCondition.create('Dirichlet', [1,2,3], 'Disp', 0)])
>>> pb.bc.add(my_bc)
For non linear problem, we can also defined the time_func to sepcify how the boundary condition evolve during time. By default, a linear evolution is enforced. The time function depend on the time_factor. The time_factor is 0 at the begining and 1 at the end of the iteration. The time function is also a factor to the prescribed value and should be between 0 and 1.
>>> def step_function(t_fact):
>>> if t_fact == 0: return 0
>>> else: return 1
>>> pb.bc.add('Dirichlet', 'nodeset', 'Disp', 1, time_func = step_function)])
Boundary conditions sharing some dof
A displacement and a generalized force may be prescribed on the same node as long as the dof are different. When several boundary conditions are defined on the same dof:
The generalized forces of the Neumann boundary conditions are summed. This includes the distributed loads (e.g.
fedoo.constraint.Pressure), whose equivalent nodal forces are computed over the whole surface as if no dof was blocked.A Dirichlet boundary condition prevails over a Neumann one: the prescribed force is ignored and the nodal force is the reaction.
For several Dirichlet boundary conditions, the last one in the list is applied.
The boundary conditions that can’t be applied together are reported the
first time they are applied to the problem: Dirichlet boundary conditions
with different values on a dof, Dirichlet boundary condition on a dof
eliminated by a multi point constraint, or dof eliminated by several multi
point constraints. A non zero nodal force ignored on a blocked dof is also
reported. By default a warning is emitted. The attribute check_mode
allows to raise an error for the incompatible boundary conditions, or to
disable the check:
>>> pb.bc.check_mode = "raise" # or "warn" (default), "ignore"
The report is also available as a list with
fedoo.ListBC.check():
>>> for record in pb.check_boundary_conditions():
>>> print(record["kind"], record["node"], record["variable"])
Multi Point Constraints
|
Class that define multi-point constraints |
Advanced BC and constraints
Fedoo also provides higher-level boundary-condition objects for common
distributed loads and constraint patterns. Distributed loads are converted to
external Neumann forces, while the kinematic constraints below are based on
MPC equations and can be added directly to pb.bc.
Distributed loads
Distributed loads are defined by dedicated assembly objects. The recommended
way to apply them is to add the load assembly to the problem boundary
conditions. If the assembly provides an as_neumann() method,
fedoo.ListBC.add() automatically converts it to an equivalent
Neumann boundary condition. This is especially useful for nonlinear problems:
the load is included in the external force vector, follows the usual load
coefficient, and is taken into account by the residual normalization.
solid_assembly = fd.Assembly.create(wf, mesh)
pressure = fd.constraint.Pressure(surface_mesh, value, nlgeom=True)
pb = fd.problem.NonLinear(solid_assembly, nlgeom=True)
pb.bc.add(pressure)
The explicit form is also available:
pb.bc.add(pressure.as_neumann())
The load assembly can still be combined directly with the mechanical assembly. This form is kept as an alternative, but is not recommended for nonlinear analyses because the load is part of the residual assembly rather than a boundary-condition force.
pb = fd.problem.Linear(solid_assembly + pressure)
|
Distributed force (e.g gravity load). |
|
Pressure load. |
|
Surface stress with a fixed orientation. |
Generic Lagrange-multiplier constraints
Any boundary-condition object that generates only fedoo.MPC leaves
can be enforced with Lagrange multipliers by wrapping it in
fedoo.LagrangeMultiplierAssembly. Original MPC
coefficients are preserved; the first term is not normalized as a slave DOF.
periodic_lm = fd.LagrangeMultiplierAssembly(
mesh,
fd.constraint.PeriodicBC(periodicity_type="small_strain"),
name="PeriodicLM",
)
pb = fd.problem.Linear(solid_assembly + periodic_lm)
Do not also add the wrapped constraint to pb.bc: that would enforce its
equations a second time through MPC elimination. The bordered system is
an indefinite saddle-point system. Use a direct solver; positive-definite
iterative solvers such as cg are not suitable.
|
Enforce MPC-generated linear constraints with Lagrange multipliers. |
Mean-value constraint
The mean value of a field can be constrained with
fedoo.constraint.MeanValueConstraint. It is an assembly object
that is summed with the stiffness assembly. Its main use is to remove rigid
body translation in periodic homogenization without pinning an arbitrary node:
pb = fd.problem.Linear(
assembly + fd.constraint.MeanValueConstraint(mesh)
)
The constraint uses one Lagrange multiplier per variable. The zero-mean condition provides a canonical periodic fluctuation field while retaining the same stress field as a single-node pin.
|
Enforce the weighted mean of a field over a set of nodes. |
Kinematic MPC constraints
Rigid tie
The fedoo.constraint.RigidTie constraint couples a set of nodes to
a rigid-body motion. It creates global degrees of freedom for rigid
translation and rotation, then eliminates the selected nodal displacement DOFs
through MPC equations. Prescribing RigidDisp or RigidRot therefore
drives the whole selected node set as a rigid object.
This constraint is useful to impose a strict rigid connection, for instance to
drive a face from a reference rigid point. It is more restrictive than
fedoo.constraint.MeanMotion, because all selected nodes follow the
same rigid-body kinematics and local warping of the selected surface is
suppressed.
|
Constraint that eliminate dof assuming a rigid body tie between nodes. |
|
Constraint that eliminate dof assuming a rigid body tie between nodes in 2D. |
Mean motion control
The fedoo.constraint.MeanMotion constraint is a less
restrictive alternative to fedoo.constraint.RigidTie. It extracts
selected components of the best-fit motion of a node set or surface mesh in a
mean sense, without forcing every selected node to follow a rigid motion. The
selected surface can therefore deform locally while selected mean translation
and rotation components are controlled.
The components argument is required and defines which global DOFs are
created. Accepted component names include "RotX", "MeanRotX",
"DispZ" and "MeanDispZ". The vector aliases "Rot"/"MeanRot"
and "Disp"/"MeanDisp" select all rotation or displacement components
available in the current modeling dimension. In 2D, the rotation vector
contains only MeanRotZ.
When a surface mesh is passed, nodal weights are computed from the element area, or from the element length in 2D. These weights define an area-weighted least-squares projection of the displacement field onto rigid-body modes.
For displacement-only mean control, select only displacement components:
mean_disp = fd.constraint.MeanMotion(surface_mesh, components="DispZ")
pb.bc.add(mean_disp)
pb.bc.add("Dirichlet", "MeanDispZ", -delta)
When finite_rotation is left to None (the default), rotational
components use finite rotations automatically if geometrical nonlinearity is
enabled on the problem. Otherwise, the rotation components are small-rotation
vector components. The finite-rotation mode uses a rotation-vector convention
and linearizes the constraint at each Newton iteration using the derivative of
the rotation matrix with respect to the rotation vector.
face_motion = fd.constraint.MeanMotion(
surface_mesh,
components="RotZ",
)
pb.bc.add(face_motion)
pb.bc.add("Dirichlet", face_motion.node_by_variable["MeanRotZ"], "MeanRotZ", angle)
face_motion = fd.constraint.MeanMotion(
surface_mesh,
components=["RotZ", "DispZ"],
)
pb.bc.add(face_motion)
pb.bc.add("Dirichlet", "MeanRotZ", angle)
pb.bc.add("Dirichlet", "MeanDispZ", -delta)
The reactions on the created global DOFs are generalized forces. For instance,
the reaction associated with MeanRotZ is the conjugate moment around the
mean-motion center.
|
Constraint defining selected global DOFs equal to a mean motion. |
Periodic boundary conditions
The fedoo.constraint.PeriodicBC constraint creates MPC equations
between opposite boundaries so that their displacement fluctuations are
periodic. It is mainly used for representative volume elements and
homogenization problems, where the boundary kinematics must be compatible with
a prescribed macroscopic strain or displacement gradient.
The constraint can build the required node pairings from opposite node sets
and then be added to pb.bc like the other MPC-based constraints.
|
Periodic boundary conditions constraint. |
Contact
Fedoo provides two contact approaches: a penalty-based method and an
IPC (Incremental Potential Contact) method. Both are implemented as
assembly objects and can be combined with other assemblies using
fedoo.Assembly.sum().
Penalty-based contact
The penalty method uses a node-to-surface formulation. It is available for 2D problems and supports frictionless contact.
|
Contact Assembly based on a node 2 surface formulation |
|
Self contact Assembly (ie contact of a geomtry between itself) based on a node 2 surface formulation |
The contact class is derived from assembly. To add a contact contraint
to a problem, we need first to create the contact assembly (using the class
fedoo.constraint.Contact) and then to add it to the global
assembly with fedoo.Assembly.sum().
IPC contact
The IPC (Incremental Potential Contact) method uses barrier potentials from the ipctk library to guarantee intersection-free configurations. It supports both 2D and 3D problems, friction, and optional CCD (Continuous Collision Detection) line search.
Unlike the penalty method, IPC does not require tuning a penalty
parameter. The barrier stiffness \(\kappa\) is automatically computed
and adaptively updated to balance the elastic and contact forces. dhat
sets the barrier activation range (default: 0.1% of the bounding-box
diagonal); dmin optionally sets a physical minimum separation.
Choosing dhat — The default dhat=1e-3 (relative) means the
barrier activates when surfaces are within 0.1 % of the bounding-box
diagonal. For problems with a very small initial gap, increase dhat
(e.g. 1e-2) so the barrier catches contact early. For tight-fitting
assemblies where a visible gap is unacceptable, decrease it (e.g.
1e-4), but expect more Newton–Raphson iterations.
Contact thickness / offset — dmin=0 preserves zero-thickness
contact. A positive dmin sets the minimum distance between collision
primitives: the barrier activates when the distance is below
dmin + dhat and diverges as it approaches dmin. dmin is always
in absolute model length units, even with dhat_is_relative=True.
Increasing dhat alone does not impose a physical thickness.
For 2D beam centerlines with a common radius of 0.2 mm, use a minimum centerline separation of 0.4 mm:
contact = fd.constraint.IPCContact(
mesh, surface_mesh=mesh,
dmin=0.4, dhat=0.05, dhat_is_relative=False,
use_ccd=True,
)
assembly = fd.Assembly.sum(beam_assembly, contact)
The initial geometry must have separation strictly greater than dmin
for every eligible collision pair. IPC excludes primitives sharing a
vertex, but nearby segments around connected joints can still require
additional collision exclusions: a vertex and the edges two rings away on
the same surface are about one edge length apart, so dmin + dhat
must stay below the edge length of the surface mesh, unless
excluded_rings is raised: the nearest primitive pair allowed to collide is
then about excluded_rings - 1 edge lengths away, which lets the shell
thickness exceed the element size.
The offset is uniform across the mesh; it does not model section deformation
or automatically add periodic-image contact. IPCSelfContact accepts the same dmin argument.
Nonzero offsets are not supported with use_ogc=True.
Exact or projected contact hessian — by default each collision hessian
is projected onto the positive semi-definite cone (psd_projection=True),
which guarantees a descent direction for the energy line search but makes
the Newton tangent inexact: the iterations then converge linearly, about 0.7
per iteration, as soon as the contact carries a large share of the load.
psd_projection=False assembles the exact hessian; with a direct solver
this restores the quadratic convergence of Newton-Raphson.
Rigid-body contact also accepts dmin through
RigidBody.set_static_obstacle(..., dmin=...) and
RigidBody.set_static_plane(..., dmin=...). All analytic planes attached
to one body must use the same dmin, dhat and kappa.
CCD line search — Enabling use_ccd=True is recommended for
problems where first contact occurs suddenly (e.g. a punch hitting a
plate) or where self-contact can cause rapid topology changes.
CCD uses the same dmin as the barrier so Newton corrections preserve
the minimum separation. Disabling CCD removes this swept-path protection.
Energy-based backtracking — When use_ccd=True, an
energy-based backtracking phase is automatically enabled after CCD
filtering: the step is halved until total energy (exact barrier +
quadratic elastic approximation) decreases. This matches the
reference IPC algorithm and improves convergence robustness. Set
line_search_energy=False to disable (faster but may degrade
convergence for difficult contact scenarios).
Convergence criterion — The 'Force' convergence criterion
is recommended for IPC contact problems. It measures the relative
decrease of the force residual, matching the gradient-norm convergence
used by reference IPC implementations:
pb.set_nr_criterion('Force', tol=5e-3, max_subiter=15)
The 'Displacement' criterion may become unreliable as contact
stiffness grows.
The ipctk package is required and can be installed with:
pip install ipctk
# or
pip install fedoo[ipc]
|
Contact Assembly based on the IPC (Incremental Potential Contact) method. |
|
Self-contact Assembly using the IPC method. |
Penalty contact example
Here an example of a contact between a square and a disk using the penalty method.
import fedoo as fd
import numpy as np
fd.ModelingSpace("2D")
filename = 'disk_rectangle_contact' #file to save results
#---- Create geometries --------------
mesh_rect = fd.mesh.rectangle_mesh(nx=11, ny=21,
x_min=0, x_max=1, y_min=0, y_max=1,
elm_type = 'quad4', name = 'Domain'
)
mesh_rect.element_sets['rect'] = np.arange(0, mesh_rect.n_elements)
mesh_disk = fd.mesh.disk_mesh(radius=0.5, nr=6, nt=6, elm_type = 'quad4')
mesh_disk.nodes+=np.array([1.5,0.48]) # translate disk on the right
mesh_disk.element_sets['disk'] = np.arange(0,mesh_disk.n_elements)
mesh = fd.Mesh.stack(mesh_rect,mesh_disk)
#node sets for boundary conditions
nodes_left = mesh.find_nodes('X',0)
nodes_right = mesh.find_nodes('X',1)
nodes_bc = mesh.find_nodes('X>1.5')
nodes_bc = list(set(nodes_bc).intersection(mesh.node_sets['boundary']))
#---- Define contact --------------
#slave surface = right face of rectangle mesh
nodes_contact = nodes_right
surf = fd.mesh.extract_surface(mesh.extract_elements('disk'))
contact = fd.constraint.Contact(nodes_contact, surf)
contact.contact_search_once = True
contact.eps_n = 5e5
#---- Material properties --------------
props = np.array([200e3, 0.3, 1e-5, 300, 1000, 0.3])
# E, nu, alpha (non used), Re, k, m
material_rect = fd.constitutivelaw.Simcoon("EPICP", props)
material_disk = fd.constitutivelaw.ElasticIsotrop(50e3, 0.3) #E, nu
material = fd.constitutivelaw.Heterogeneous(
(material_rect, material_disk),
('rect', 'disk')
)
#---- Build problem --------------
wf = fd.weakform.StressEquilibrium(material, nlgeom = True)
solid_assembly = fd.Assembly.create(wf, mesh)
assembly = fd.Assembly.sum(solid_assembly, contact)
pb = fd.problem.NonLinear(assembly)
results = pb.add_output(filename,
solid_assembly,
['Disp', 'Stress', 'Strain', 'Fext']
)
pb.bc.add('Dirichlet',nodes_left, 'Disp',0)
pb.bc.add('Dirichlet',nodes_bc, 'Disp', [-0.4,0.2])
pb.set_nr_criterion("Displacement", tol = 5e-3, max_subiter = 5)
#---- Solve problem in two steps: load, unload --------------
pb.nlsolve(dt = 0.005, tmax = 1, update_dt = True, interval_output = 0.01)
pb.bc.remove(-1) #remove last boundary contidion
pb.bc.add('Dirichlet',nodes_bc, 'Disp', [0,0])
pb.nlsolve(dt = 0.005, tmax = 1, update_dt = True, interval_output = 0.01)
# =============================================================
# Example of plots with pyvista - uncomment the desired plot
# =============================================================
# ------------------------------------
# Simple plot with default options
# ------------------------------------
results.plot('Stress', component='vm', data_type='Node')
results.plot('Disp', component = 0, data_type='Node')
# ------------------------------------
# Write movie with default options
# ------------------------------------
results.write_movie(filename,
'Stress',
component = 'XX',
data_type = 'Node',
framerate = 24,
quality = 5,
clim = [-3e3, 3e3]
)
Video of results:
contact video
IPC contact example
The same disk-rectangle contact problem can be solved with the IPC method. The IPC method does not require choosing slave/master nodes or tuning a penalty parameter. The barrier stiffness is automatically computed and adapted.
import fedoo as fd
import numpy as np
fd.ModelingSpace("2D")
#---- Create geometries (same as penalty example) --------------
mesh_rect = fd.mesh.rectangle_mesh(nx=11, ny=21,
x_min=0, x_max=1, y_min=0, y_max=1,
elm_type='quad4', name='Domain')
mesh_disk = fd.mesh.disk_mesh(radius=0.5, nr=6, nt=6, elm_type='quad4')
mesh_disk.nodes += np.array([1.5, 0.48])
mesh = fd.Mesh.stack(mesh_rect, mesh_disk)
#---- Define IPC contact --------------
surf = fd.mesh.extract_surface(mesh)
ipc_contact = fd.constraint.IPCContact(
mesh, surface_mesh=surf,
friction_coefficient=0.3, # Coulomb friction coefficient
use_ccd=True, # enable CCD line search for robustness
)
# barrier_stiffness is auto-computed; dhat defaults to 1e-3 * bbox_diag
#---- Material and problem setup --------------
material = fd.constitutivelaw.ElasticIsotrop(200e3, 0.3)
wf = fd.weakform.StressEquilibrium(material, nlgeom=True)
solid_assembly = fd.Assembly.create(wf, mesh)
assembly = fd.Assembly.sum(solid_assembly, ipc_contact)
pb = fd.problem.NonLinear(assembly)
res = pb.add_output('results', solid_assembly, ['Disp', 'Stress'])
# ... add BCs ...
pb.nlsolve(dt=0.005, tmax=1)
Note
When using add_output, pass the solid assembly (not the sum).
AssemblySum objects cannot be used directly for output extraction.
IPC self-contact example
For self-contact problems, use IPCSelfContact
which automatically extracts the surface from the volumetric mesh.
import fedoo as fd
import numpy as np
fd.ModelingSpace("3D")
mesh = fd.Mesh.read("gyroid.vtk")
material = fd.constitutivelaw.ElasticIsotrop(1e5, 0.3)
# Self-contact: auto surface extraction, auto barrier stiffness
# Add line_search_energy=True for extra robustness (slower)
contact = fd.constraint.IPCSelfContact(mesh, use_ccd=True)
wf = fd.weakform.StressEquilibrium(material, nlgeom="UL")
solid = fd.Assembly.create(wf, mesh)
assembly = fd.Assembly.sum(solid, contact)
pb = fd.problem.NonLinear(assembly)
res = pb.add_output("results", solid, ["Disp", "Stress"])
# ... add BCs ...
pb.nlsolve(dt=0.05, tmax=1, update_dt=True)