Constraints
Constraints are defined globally, and may be used for introducing additional scalar equations into the global equation system.
Example syntax:
*constraint, type=rigidbody, name=myRigigBody
nSet=right
referencePoint=rBottom
equalvaluelagrangian - Constrain nodal values to equal values
Module edelweissfe.constraints.equalvaluelagrangian
Option |
Type / Default |
Description |
|---|---|---|
field |
str (required) |
The field this constraint acts on. |
component |
int (required) |
The component of the field. |
nSet |
str (required) |
The node set to be constrained. |
testfiles/marmot/EqualValueLagrangianConstraint/test.inp*modelGenerator, generator=planeRectQuad, name=gen
elType=CPS8R
nX=20
nY=20
l=100
h=100
*constraint, type=equalValueLagrangian, name=myConstraint
nSet=gen_top
field=displacement
component=1
*material, name=linearelastic, id=linearelastic
30000.0, 0.15
*section, name=section1, thickness=1.0, material=linearelastic, type=plane
all
*job, name=myJob, domain=2d
*solver, solver=NIST, name=theSolver
*fieldOutput
>>perNode, name=displacement, elSet=all, field=displacement, result=U,
>>perNode, name=RF, nSet=gen_bottom, field=displacement, result=P,
>>perElement, name=stress, elSet=all, result=stress, quadraturePoint=0:4, f(x)='mean(x, axis=1)'
>>perElement, name=strain, elSet=all, result=strain, quadraturePoint=0:4, f(x)='mean(x, axis=1)'
*output, type=monitor, name=theMonitor
fieldOutput = RF, f(x) = 'sum ( x, axis=0)'
*output, type=ensight, name=ensightExport
>>perNode, fieldOutput=displacement
>>perElement, fieldOutput=stress
>>perElement, fieldOutput=strain
>>configuration, overwrite=yes
*step, solver=theSolver
maxInc=1e-1, minInc=1e-5, maxNumInc=1000, maxIter=25, stepLength=1
>>options, name=theSolver, extrapolation=linear
>>dirichlet, name=bottom, nSet=gen_bottom, field=displacement, 1=0.0, 2=0.0
>>dirichlet, name=load, nSet=gen_rightTop, field=displacement, 1=5.0, 2=5.0
equalvaluepenalty - Constrain nodal values to equal values
Module edelweissfe.constraints.equalvaluepenalty
Option |
Type / Default |
Description |
|---|---|---|
field |
str (required) |
The field this constraint acts on. |
component |
int (required) |
The component of the field. |
penalty |
float (required) |
The numerical penalty value. |
nSet |
str (required) |
The node set to be constrained. |
testfiles/marmot/EqualValuePenaltyConstraint/test.inp*modelGenerator, generator=planeRectQuad, name=gen
elType=CPS8R
nX=20
nY=20
l=100
h=100
*constraint, type=equalValuePenalty, name=myConstraint
nSet=gen_top
field=displacement
component=1
penalty=1e8
*material, name=linearelastic, id=linearelastic
30000.0, 0.15
*section, name=section1, thickness=1.0, material=linearelastic, type=plane
all
*job, name=myJob, domain=2d
*solver, solver=NIST, name=theSolver
*fieldOutput
>>perNode, name=displacement, elSet=all, field=displacement, result=U,
>>perNode, name=RF, nSet=gen_bottom, field=displacement, result=P,
>>perElement, name=stress, elSet=all, result=stress, quadraturePoint=0:4, f(x)='mean(x, axis=1)'
>>perElement, name=strain, elSet=all, result=strain, quadraturePoint=0:4, f(x)='mean(x, axis=1)'
*output, type=monitor, name=theMonitor
fieldOutput = RF, f(x) = 'sum ( x, axis=0)'
*output, type=ensight, name=ensightExport
>>perNode, fieldOutput=displacement
>>perElement, fieldOutput=stress
>>perElement, fieldOutput=strain
>>configuration, overwrite=yes
*step, solver=theSolver
maxInc=1e-1, minInc=1e-5, maxNumInc=1000, maxIter=25, stepLength=1
>>options, name=theSolver, extrapolation=linear
>>dirichlet, name=bottom, nSet=gen_bottom, field=displacement, 1=0.0, 2=0.0
>>dirichlet, name=load, nSet=gen_rightTop, field=displacement, 1=5.0, 2=5.0
linearizedrigidbody - Linearized rigid body constraints
Module edelweissfe.constraints.linearizedrigidbody
Option |
Type / Default |
Description |
|---|---|---|
nSet |
str (required) |
Node set to tie. |
referencePoint |
str (required) |
Node set containing only the reference point. |
testfiles/marmot/LinearizedRigidBodyConstraint/test.inp*modelGenerator, generator=planeRectQuad, name=gen
l=100
h=10
nX=15
nY=5
elType=CPS4
*material, name=linearElastic, id=le
30000.0, 0.15
*section, name=section1, thickness=1.0, material=le, type=plane
all
*job, name=cps4job, domain=2d
*solver, solver=NIST, name=theSolver
*fieldOutput
>>perNode, name=load, nSet=gen_left, field=displacement, result=P, f(x)='sum(x[:,1])'
>>perNode, name=displacement, elSet=all, field=displacement, result=U
*output, type=monitor
fieldOutput=load
*constraint, type=linearizedRigidBody, name=rb1
nSet=gen_right
referencePoint=gen_rightTop
*output, type=ensight, name=ensightExport
>>perNode, fieldOutput=displacement
>>configuration, overwrite=yes
*step, solver=theSolver
maxInc=1e-1, minInc=1e-4, maxNumInc=1000, maxIter=4
>>options, name=theSolver, extrapolation=linear
>>dirichlet, name=leftCorner, nSet=gen_leftBottom, field=displacement, 1=0.0, 2=0.0
>>dirichlet, name=rotRP, nSet=gen_rightTop, field=rotation, 1=0.05
>>dirichlet, name=movRP, nSet=gen_rightTop, field=displacement, 2=0.00
testfiles/marmot/LinearizedRigidBodyConstraint2D/test.inp*modelGenerator, generator=planeRectQuad, name=gen
l=100
h=10
nX=15
nY=5
elType=CPS4
*node, nSet=RP
999999, 150, 5
*material, name=linearElastic, id=myMaterial
30000.0, 0.15
*section, name=section1, thickness=1.0, material=myMaterial, type=plane
all
*job, name=myJob, domain=2d
*solver, solver=NIST, name=theSolver
*fieldOutput
>>perNode, name=load, nSet=gen_left, field=displacement, result=P, f(x)='sum(x[:,1])'
>>perNode, name=displacement, elSet=all, field=displacement, result=U
*output, type=monitor
fieldOutput=load
*constraint, type=linearizedRigidBody, name=fixedSupport
nSet=gen_left
referencePoint=gen_leftTop
*constraint, type=linearizedRigidBody, name=pinnedSupport
nSet=gen_right
referencePoint=RP
*output, type=ensight, name=test
>>perNode, fieldOutput=displacement
>>configuration, overwrite=yes
*step, solver=theSolver
maxInc=1e-1, minInc=1e-4, maxNumInc=1000, maxIter=4
>>options, name=theSolver, extrapolation=linear
>>dirichlet, name=rotFixedSupport, nSet=gen_leftTop, field=rotation, 1=0
>>dirichlet, name=movFixedSupport, nSet=gen_leftTop, field=displacement, 1=0, 2=0
>>dirichlet, name=movPinnedSupport, nSet=RP, field=displacement, 1=0, 2=0
*step, solver=theSolver
maxInc=1e-1, minInc=1e-4, maxNumInc=1000, maxIter=4
>>options, name=theSolver, extrapolation=linear
>>distributedload, name=loadTop, surface=gen_top, type=pressure, magnitude=10.0, f(t)=1
testfiles/marmot/LinearizedRigidBodyConstraint3D/test.inp*modelGenerator, generator=boxGen, name=gen
lX=100
lY=10
lZ=10
nX=15
nY=5
nZ=5
elType=C3D8
*node, nSet=RP
999999, 150, 5, 5
*material, name=linearElastic, id=myMaterial
30000.0, 0.15
*section, name=section1, material=myMaterial, type=solid
all
*job, name=myJob, domain=3d
*solver, solver=NIST, name=theSolver
*fieldOutput
>>perNode, name=load, nSet=gen_left, field=displacement, result=P, f(x)='sum(x[:,1])'
>>perNode, name=displacement, elSet=all, field=displacement, result=U
*output, type=monitor
fieldOutput=load
*constraint, type=linearizedRigidBody, name=fixedSupport
nSet=gen_left
referencePoint=gen_topLeftBack
*constraint, type=linearizedRigidBody, name=pinnedSupport
nSet=gen_right
referencePoint=RP
*output, type=ensight, name=test
>>perNode, fieldOutput=displacement
>>configuration, overwrite=yes
*step, solver=theSolver
maxInc=1e-1, minInc=1e-4, maxNumInc=1000, maxIter=4
>>options, name=theSolver, extrapolation=linear
>>dirichlet, name=rotFixedSupport, nSet=gen_topLeftBack, field=rotation, 1=0, 2=0, 3=0
>>dirichlet, name=movFixedSupport, nSet=gen_topLeftBack, field=displacement, 1=0, 2=0, 3=0
>>dirichlet, name=movPinnedSupport, nSet=RP, field=displacement, 1=0, 2=0, 3=0
*step, solver=theSolver
maxInc=1e-1, minInc=1e-4, maxNumInc=1000, maxIter=4
>>options, name=theSolver, extrapolation=linear
>>distributedload, name=loadTop, surface=gen_top, type=pressure, magnitude=10.0, f(t)=1
>>dirichlet, name=rotation, nSet=RP, field=rotation, 2=0.1
linearizedrigidbodympc - Linearized rigid body constraints by DOF elimination
Module edelweissfe.constraints.linearizedrigidbodympc
Linearized rigid body constraint enforced by master-slave DOF elimination (multi-point constraint):
the kinematics of linearizedrigidbody, without its Lagrange
multipliers, so without extra DOFs and without a saddle-point system (iterative linear solvers such
as blockamg need a nonzero diagonal).
The displacement of every slave node \(s\) follows the displacement and the (small) rotation of the reference point (RP),
in 3D (rotation field with 3 components), and in 2D (one rotation \(\phi\) about z)
The slave displacements are eliminated; the displacement and the rotation of the RP are the masters. Prescribe Dirichlet conditions on the RP only, not on the slave nodes: a slave degree of freedom that is also prescribed loses its record to the Dirichlet condition (reported by the solver), so the body is no longer rigid at that node.
This constraint is a MeshDependent: if adaptive refinement
changes the slave node set (e.g. by hadaptivity), the
records are rebuilt from the current members, so new nodes are tied as well.
Limitation: static solvers only. The RP carries no mass of its own, and both dynamic solvers check
the assembled inertia before the multi-point constraint enters, so they reject the RP’s degrees of
freedom (NID: no mass on a time-integrated DOF, as for
linearizedrigidbody; NED: an inertia on only some DOFs of a field).
For NED, folding the slave masses onto the RP would not suffice: the row-sum lumping used for ties
gives the RP rotation the signed static moment of the slave masses instead of a rotational inertia.
Option |
Type / Default |
Description |
|---|---|---|
nSet |
str (required) |
Node set to tie. |
referencePoint |
str (required) |
Node set containing only the reference point. |
testfiles/marmot/LinearizedRigidBodyMPCConstraint3D/test.inp*modelGenerator, generator=boxGen, name=gen
lX=100
lY=10
lZ=10
nX=15
nY=5
nZ=5
elType=C3D8
*node, nSet=RP
999999, 150, 5, 5
*material, name=linearElastic, id=myMaterial
30000.0, 0.15
*section, name=section1, material=myMaterial, type=solid
all
*job, name=myJob, domain=3d
*solver, solver=NIST, name=theSolver
*fieldOutput
>>perNode, name=load, nSet=gen_left, field=displacement, result=P, f(x)='sum(x[:,1])'
>>perNode, name=displacement, elSet=all, field=displacement, result=U
*output, type=monitor
fieldOutput=load
*constraint, type=linearizedRigidBodyMPC, name=fixedSupport
nSet=gen_left
referencePoint=gen_topLeftBack
*constraint, type=linearizedRigidBodyMPC, name=pinnedSupport
nSet=gen_right
referencePoint=RP
*output, type=ensight, name=test
>>perNode, fieldOutput=displacement
>>configuration, overwrite=yes
*step, solver=theSolver
maxInc=1e-1, minInc=1e-4, maxNumInc=1000, maxIter=4
>>options, name=theSolver, extrapolation=linear
>>dirichlet, name=rotFixedSupport, nSet=gen_topLeftBack, field=rotation, 1=0, 2=0, 3=0
>>dirichlet, name=movFixedSupport, nSet=gen_topLeftBack, field=displacement, 1=0, 2=0, 3=0
>>dirichlet, name=movPinnedSupport, nSet=RP, field=displacement, 1=0, 2=0, 3=0
*step, solver=theSolver
maxInc=1e-1, minInc=1e-4, maxNumInc=1000, maxIter=4
>>options, name=theSolver, extrapolation=linear
>>distributedload, name=loadTop, surface=gen_top, type=pressure, magnitude=10.0, f(t)=1
>>dirichlet, name=rotation, nSet=RP, field=rotation, 2=0.1
rigidbody - Geometrically exact rigid body constraints in 3D
Module edelweissfe.constraints.rigidbody
Option |
Type / Default |
Description |
|---|---|---|
nSet |
str (required) |
Node set to tie. |
referencePoint |
str (required) |
Node set containing only the reference point. |
testfiles/marmot/RigidBodyConstraintLargeDeformations3D/test.inp*include, input=mesh.inp
*material, name=COMPRESSIBLENEOHOOKE, id=le,
3500, 1500
*section, name=section1, material=le, type=solid
all
*job, name=c3d8job, domain=3d
*solver, solver=NISTParallel, name=theSolver
*fieldOutput
>>perNode, name=displacement, elSet=all, field=displacement, result=U
>>perElement, name=stress, elSet=all, result=stress, f(x)=' mean(x, axis=1) ', quadraturePoint=0:8
*constraint, type=rigidbody, name=rb1
nSet=right
referencePoint=rBottom
*output, type=ensight, name=ensightExport
>>perNode, fieldOutput=displacement
>>perElement, fieldOutput=stress
>>configuration, overwrite=yes
*step, solver=theSolver
maxInc=5e-2, minInc=1e-4, maxNumInc=1000, maxIter=10
>>options, name=theSolver, extrapolation=linear
>>dirichlet, name=left, nSet=left, field=displacement, 1=0.0, 2=0, 3=0
>>dirichlet, name=rotRP, nSet=rBottom, field=rotation, 1=3.141,
penaltyindirectcontrol - Penalty based indirect control
Module edelweissfe.constraints.penaltyindirectcontrol
Option |
Type / Default |
Description |
|---|---|---|
field |
str, default=’displacement’ |
The field this constraint acts on. |
cVector |
str (required) |
The projection vector for the constrained nodes (e.g., CMOD). |
constrainedNSet |
str (required) |
The node set for determining the constraint (e.g., CMOD). |
loadNSet |
str (required) |
The node set for application of the controlled load. |
loadVector |
str (required) |
The vector (in correct) dimensions and tensorial order determining the load. |
length |
float (required) |
The value of the constraint (e.g., CMOD). |
penaltyStiffness |
float (required) |
The stiffness for formulating the constraint. |
offset |
float, default=0.0 |
A correction value for the computation of the constraint (e.g, initial displacement). |
normalizeLoad |
bool, default=True |
Normalize the applied force per node w. r. t. the number of nodes, i.e., apply a load irrespective of the total number of nodes in ``loadNSet``. |
f(t) |
str, default=None |
Amplitude function. |
testfiles/marmot/PenaltyBasedIndirectControl/test.inp*material, name=LinearElastic, id=myMaterial
**Isotropic
**E nu
1000, 0.2
*job, name=job, domain=2d
*solver, solver=NIST, name=theSolver
*modelGenerator, generator=planerectquad, name=gen
nX =5
nY =5
l =20
h =50
elType = CPE4
*section, name=section1, material=myMaterial, type=plane, thickness=1
all
*fieldOutput
>>perNode, elSet=all, field=displacement, result=U, name=displacement
>>perNode, nSet=gen_top, field=displacement, result=U, name=disp_top
>>perNode, nSet=gen_bottom, field=displacement, result=P, name=RF
>>perElement, elSet=all, result=strain, name=strain, quadraturePoint=0:4, f(x)='np.mean(x,axis=1)'
>>perElement, elSet=all, result=stress, name=stress, quadraturePoint=0:4, f(x)='np.mean(x,axis=1)'
*output, type=ensight,
>>perNode, fieldOutput=displacement
>>perElement, fieldOutput=stress
>>perElement, fieldOutput=strain
>>configuration, overwrite=yes
*output, type=monitor, name=dispMon
fieldOutput=disp_top, f(x)='(x[:,1])'
fieldOutput=RF, f(x)='sum(x[:,1])'
*nset, nset=measurement_nodes
gen_leftBottom
gen_leftTop
*constraint, type=penaltyindirectcontrol, name=pc1
field=displacement
cVector='0, 1, 0, -1'
constrainedNSet=measurement_nodes
loadNSet=gen_rightTop
loadVector='0, -1'
length=-1
penaltyStiffness=1e8
offset=0
*step, solver=theSolver
maxInc=1e0, minInc=1e0, maxNumInc=100, maxIter=25, stepLength=1
>>modelupdate, update='model.constraints["pc1"].active=False'
>>dirichlet, name=fix_bottom, nSet=gen_bottom, field=displacement, 1=0, 2=0
*step, solver=theSolver
maxInc=2e-2, minInc=1e-6, maxNumInc=1000, maxIter=25, stepLength=1
>>modelupdate, update='model.constraints["pc1"].active=True'
directionalspringpenalty - Assigning a stiffness to specific degrees of freedom
Module edelweissfe.constraints.directionalspringpenalty
Option |
Type / Default |
Description |
|---|---|---|
field |
str (required) |
The field this constraint acts on. |
component |
int (required) |
The component of the field. |
penalty |
float (required) |
The numerical penalty value. |
nSet |
str (required) |
The node set to be constrained. |
testfiles/marmot/DirectionalSpringPenaltyConstraint/test.inp*modelGenerator, generator=boxGen, name=gen
nX =1
nY =1
nZ =1
lX =1
lY =1
lZ =1
elType =C3D8
*material, name=LinearElastic, id=myMaterial
3000000, 0.2
*section, name=section1, material=myMaterial, type=solid
all
*job, name=myJob, domain=3d
*solver, name=theSolver, solver=NISTParallel
*fieldOutput
** for ensight output
>>perNode, elSet=all, name=displacement, field=displacement, result=U
>>perElement, elSet=all, name=stress, result=stress, quadraturePoint=0:8, f(x)='mean(x, axis=1)'
>>perElement, elSet=all, name=strain, result=strain, quadraturePoint=0:8, f(x)='mean(x, axis=1)'
** for xy-data
>>perNode, name=RF, nSet=gen_right, field=displacement, result=P, f(x)='-sum(x[:, 0])'
>>perNode, name=Uleft, nSet=gen_left, field=displacement, result=U, f(x)='mean(x[:, 0])'
*output, type=monitor, name=RFMon
fieldOutput=RF
fieldOutput=Uleft
*output, type=ensight, name=ensightExport
>>perNode, fieldOutput=displacement
>>perElement, fieldOutput=strain
>>perElement, fieldOutput=stress
>>configuration, overwrite=yes
*constraint, type=directionalspringpenalty, name=x_constraintLeftTop
nSet=gen_left
field=displacement
component=0
penalty=1e2
*step, solver=theSolver
maxInc=1.0, minInc=1e-5, maxNumInc=10000, maxIter=25, stepLength=1
>>dirichlet, name=leftFixed, nSet=gen_bottomLeftFront, field=displacement, 2=0., 3=0.
>>dirichlet, name=leftRot, nSet=gen_bottomLeftBack, field=displacement, 2=0.
*step, solver=theSolver
maxInc=1e-1, minInc=1e-6, maxNumInc=100000, maxIter=25, stepLength=1
>>options, name=theSolver, extrapolation=linear
>>dirichlet, name=right, nSet=gen_right, field=displacement, 1=0.1
*step, solver=theSolver
maxInc=1e-1, minInc=1e-6, maxNumInc=100000, maxIter=25, stepLength=1
>>options, name=theSolver, extrapolation=linear
>>dirichlet, name=right, nSet=gen_right, field=displacement, 1=-0.15
nodetorigidsurfacepenalty - Preventing nodes from penetrating a defined rigid boundary
Module edelweissfe.constraints.nodetorigidsurfacepenalty
Option |
Type / Default |
Description |
|---|---|---|
field |
str (required) |
The field this constraint acts on. |
component |
int (required) |
The component of the field. |
penalty |
float (required) |
The numerical penalty value. |
nSet |
str (required) |
The node set to be constrained. |
value |
float, default=0.0 |
The prescribed distance to the rigid boundary. A value of 0.0 implies no initial gap between the node set and the boundary. |
direction |
float, default=1.0 |
The normal direction outward from the continuum towards the boundary (1.0 or -1.0). |
type |
str, default=’linear’ |
The formulation type: ‘linear’ (linear force, constant stiffness with jump) or ‘quadratic’ (quadratic force, linear stiffness). |
testfiles/marmot/NodeToRigidSurfacePenaltyConstraintLinear/test.inp*modelGenerator, generator=boxGen, name=gen
nX =1
nY =1
nZ =1
lX =1
lY =1
lZ =1
elType =C3D8
*material, name=LinearElastic, id=myMaterial
3000000, 0.2
*section, name=section1, material=myMaterial, type=solid
all
*job, name=myJob, domain=3d
*solver, name=theSolver, solver=NISTParallel
*fieldOutput
** for ensight output
>>perNode, elSet=all, name=displacement, field=displacement, result=U,
>>perElement, elSet=all, name=stress, result=stress, quadraturePoint=0:1
>>perElement, elSet=all, name=strain, result=strain, quadraturePoint=0:1
** for xy-data
>>perNode, name=RF, nSet=gen_right, field=displacement, result=P, f(x)='-sum(x[:, 0])',
>>perNode, name=Uleft, nSet=gen_left, field=displacement, result=U, f(x)='mean(x[:, 0])',
*output, type=monitor, name=RFMon
fieldOutput=RF
fieldOutput=Uleft
*output, type=ensight, name=ensight_tens
>>perNode, fieldOutput=displacement
>>perElement, fieldOutput=strain
>>perElement, fieldOutput=stress
>>configuration, overwrite=yes
*constraint, type=nodetorigidsurfacepenalty, name=x_constraintLeft
nSet=gen_left
field=displacement
component=0
penalty=1e3
value=0.0
direction=-1
type=linear
*step, solver=theSolver
maxInc=1.0, minInc=1e-5, maxNumInc=10000, maxIter=25, stepLength=1
>>dirichlet, name=leftFixed, nSet=gen_bottomLeftFront, field=displacement, 2=0., 3=0.
>>dirichlet, name=leftRot, nSet=gen_bottomLeftBack, field=displacement, 2=0.
*step, solver=theSolver
maxInc=1e-1, minInc=1e-6, maxNumInc=100000, maxIter=25, stepLength=1
>>options, name=ensight_tens, intermediateSaveInterval=2, minDTForOutput=0.01
>>options, name=theSolver, extrapolation=linear
>>dirichlet, name=right, nSet=gen_right, field=displacement, 1=-0.1
*step, solver=theSolver
maxInc=1e-1, minInc=1e-6, maxNumInc=100000, maxIter=25, stepLength=1
>>options, name=ensight_tens, intermediateSaveInterval=2, minDTForOutput=0.01
>>options, name=theSolver, extrapolation=linear
>>dirichlet, name=right, nSet=gen_right, field=displacement, 1=0.2
nodetodiscreterigidbodypenalty - Contact against a discrete rigid body (deprecated)
Deprecated since version 26.11: Use surfacetodiscreterigidbodypenalty below. On quadratic (hexa20/quad8) slave faces a
node-based penalty cannot transmit the tensile corner loads of a uniform pressure, so the corner
nodes lift off (see Measured: the discrete solution lifts the corners off instead), and its per-node penalty makes the contact
stiffness depend on the mesh.
Module edelweissfe.constraints.nodetodiscreterigidbodypenalty
This is a frictionless, penalty-based, unilateral (contact-only) constraint that prevents the nodes
of a slave node set from penetrating the moving surface of a
discrete rigid body. It is the discrete-mesh counterpart of
nodetorigidsurfacepenalty (which uses an analytical boundary) and reuses the rigid body’s
own closest-point surface query for the geometry evaluation. Currently 3D only, and implicit /
explicit-static only (see Solvers).
Theory
Let a slave node \(s\) have reference position \(\mathbf{X}_s\) and displacement DOF \(\mathbf{u}_s\), so its current position is \(\mathbf{x}_s = \mathbf{X}_s + \mathbf{u}_s\). The rigid body carries a reference point (RP) with displacement DOF \(\mathbf{u}_{RP}\) and a rotation pseudo-vector DOF \(\boldsymbol{\theta}_{RP}\); its current position is \(\mathbf{x}_{RP} = \mathbf{X}_{RP} + \mathbf{u}_{RP}\) and its orientation is \(\mathbf{R} = \exp(\mathrm{skew}(\boldsymbol{\theta}_{RP}))\) (see Rigid bodies for the rigid-body kinematics).
For each slave node the rigid body’s surface query returns the signed distance \(d_s\) (negative when the node is inside the body) and the outward unit normal \(\mathbf{n}_s\) of the closest surface point, evaluated for the body’s current pose. Contact is active whenever \(d_s < 0\), with gap (penetration)
The penalty normal force is, for the two available type formulations,
with \(k\) the user-supplied penalty.
With the contact-point moment arm \(\mathbf{r}_s = \mathbf{x}_s - \mathbf{x}_{RP}\) about the current RP, the (negative) gap gradient with respect to the involved DOFs \((\mathbf{u}_s,\, \mathbf{u}_{RP},\, \boldsymbol{\theta}_{RP})\) defines the generalized direction vectors
and, since the contact force on the nodes is \(-f_n\, \mathbf{w}\), the residual contribution – its negative, see The residual and the tangent – and its tangent are assembled exactly as in the node-to-rigid-surface penalty constraint,
per slave node, where \(\mathbf{w}\) collects the slave and RP blocks.
The factor \(\mathbf{J}_r(\boldsymbol{\theta}_{RP})\) is the right Jacobian of the SO(3) exponential map,
It is required because the stored rotation DOF \(\boldsymbol{\theta}_{RP}\) is the raw total pseudo-vector (fed through the exponential map every iteration), and the exponential map does not compose additively: incrementing \(\boldsymbol{\theta}_{RP}\) by \(\delta\boldsymbol{\theta}\) does not produce a physical spin of \(\delta\boldsymbol{\theta}\) unless the body is currently unrotated. The map \(\mathbf{R}\,\mathbf{J}_r\) converts a perturbation of the stored DOF into the physical spatial rotation it actually produces; it reduces to the identity at \(\boldsymbol{\theta}_{RP} = \mathbf{0}\), where \(\mathbf{w}_{RP}\) collapses to the naive \(\mathbf{r}_s \times \mathbf{n}_s\). Omitting it makes the RP-moment term a few percent off at \(\sim 0.1\) rad accumulated rotation and \(\sim 15\%\) off at \(\sim 0.4\) rad.
Both \(\mathbf{n}_s\) and \(\mathbf{r}_s\) are recomputed from the current, total
solution every Newton iteration (no per-increment caching), so the coupling between the RP DOFs and
the contact residual is exact. Because the rigid surface is a triangulated (piecewise-planar) mesh,
\(\mathbf{n}_s\) is exactly constant within a facet — the intra-facet curvature term of the
tangent vanishes identically — but it still rotates rigidly with the RP. The geometric-stiffness
term \(f_n\, \nabla \mathbf{w}\) accounts for that rigid rotation of the normal and for the
\(\mathbf{J}_r\) dependence, giving the exact linearization of \(\mathbf{w}\) in the
implemented blocks (verified against finite differences for arbitrary accumulated rotation). The
\(\boldsymbol{\theta}_{RP}\boldsymbol{\theta}_{RP}\) self-block, which would require the second
derivative of the exponential map, is omitted; this affects only the convergence rate, not the
correctness of the converged solution. The per-block breakdown of the tangent is given in the inline
comments of the computeKernels method in the module source.
Option |
Type / Default |
Description |
|---|---|---|
nSet |
str (required) |
The (slave) node set to be protected from penetrating the rigid body. |
rigidBody |
str (required) |
The name of the discrete rigid body (as registered in model.rigidBodies). |
penalty |
float (required) |
The numerical penalty value. |
type |
str, default=’linear’ |
The formulation type: ‘linear’ (linear force, constant stiffness with jump) or ‘quadratic’ (quadratic force, linear stiffness). |
searchDistance |
float, default=None |
An optional broadphase distance (passed on to the rigid body’s surface query) for culling nodes far away from the rigid body. If not given, every slave node is queried exactly every iteration. |
testfiles/edelweiss-only/NodeToDiscreteRigidBodyContact/test.inp*material, name=LinearElastic, id=linearelastic, provider=edelweiss
**Isotropic
**E nu
1.8e4, 0.22
*job, name=c3d8job, domain=3d
*solver, solver=NIST, name=theSolver
*modelGenerator, generator=boxGen, name=block
nX =12
nY =12
nZ =2
lX =10
lY =10
lZ =5
elProvider =edelweiss
elType =C3D8
** A cylinder (radius 5, axis along z, flat end caps) taken from
** EdelweissMeshfree's discrete rigid body contact example, translated so its
** top cap sits flush with the block's bottom (z=0) face and shifted
** off-center relative to the block's footprint -- so only part of the
** cylinder's flat cap (and, at the far corner, its rounded rim) underlies
** the block, making the contact patch asymmetric.
*modelGenerator, generator=discreteRigidBodyGenerator, name=cylinder
filename=rigid_body.exo
translation='0, -10, -11.0'
** Scoped to block_all rather than all: all also contains the rigid body's
** surface/RP nodes/elements, which cannot take a solid section.
*section, name=section1, material=linearelastic, type=solid
block_all
*constraint, type=nodeToDiscreteRigidBodyPenalty, name=contact
nSet=block_back
rigidBody=cylinder
penalty=1e8
*fieldOutput
>>perNode, elSet=all, field=displacement, name=displacement, result=U
>>perNode, rigidBody=cylinder, field=displacement, name=cylinderDisplacement, result=U
>>perNode, nSet=cylinder_rp, field=displacement, name=rpForce, result=P, f(x)=' sum ( x, axis=0 )'
*output, type=ensight, name=test
>>perNode, fieldOutput=displacement
>>perNode, fieldOutput=cylinderDisplacement
>>configuration, overwrite=yes
*output, type=monitor, name=rpForceMon
fieldOutput=rpForce
*step, solver=theSolver
maxInc=0.1, minInc=1e-9, maxNumInc=100, maxIter=25, stepLength=1
>>options, name=theSolver, extrapolation=off
>>dirichlet, name=frontFixed, nSet=block_front, field=displacement, 1=0.0, 2=0.0, 3=0.0
>>dirichlet, name=rpDisp, nSet=cylinder_rp, field=displacement, 1=0.0, 2=0.0, 3=3.00
>>dirichlet, name=rpRot, nSet=cylinder_rp, field=rotation, 1=0.15, 2=0.0, 3=0.0
nodetodeformablesurfacepenalty - Node-to-deformable-surface contact
The complete theory of this constraint – surface discretization and triangulation, tributary areas, finite- and small-sliding gap kinematics, penalty and augmented-Lagrange normal laws, Coulomb friction with consistent tangents, solver integration, and the verification methodology – is documented in Contact mechanics: theory.
Module edelweissfe.constraints.nodetodeformablesurfacepenalty
Option |
Type / Default |
Description |
|---|---|---|
slaveSurface |
str (required) |
The element set of contact facet elements (Tria3ContactFacet/Line2ContactFacet) forming the slave surface; its nodes act as the contact points, penalty-weighted by their tributary areas. |
masterSurface |
str (required) |
The element set of contact facet elements (Tria3ContactFacet/Line2ContactFacet) forming the master surface. |
penalty |
float (required) |
The numerical penalty value, an interface stiffness modulus per unit slave surface area. |
type |
str, default=’linear’ |
The formulation type: ‘linear’ (linear force, constant stiffness with jump) or ‘quadratic’ (quadratic force, linear stiffness). |
searchDistance |
float, default=None |
An optional broadphase distance for the per-increment candidate-facet search. If not given, every slave is always assigned its single closest facet, without a distance gate. |
sliding |
str, default=’finite’ |
The kinematic treatment of the contact geometry: ‘finite’ (gap, gradient and exact Hessian recomputed from the current Newton iterate every iteration) or ‘small’ (Abaqus-style small sliding: the closest-point projection – facet, clamped local coordinates, and normal – is frozen once per increment from the last converged configuration, making the gap linear in the displacement DOFs). |
mu |
float, default=0.0 |
The Coulomb friction coefficient. Requires sliding=small; mu=0 disables friction. Strongly recommended in combination with type=quadratic: the quadratic law’s contact stiffness vanishes continuously at gap activation, keeping the frictional tangent continuous for slave nodes lifting off/touching down – with type=linear, the activation stiffness jump scaled by mu makes Newton prone to limit-cycling at such events. |
tangentPenalty |
float, default=None |
The tangential penalty stiffness per unit slave surface area for frictional stick. Defaults to the normal penalty. |
augmentedLagrange |
bool, default=False |
Augment the penalty force with a per-slave normal traction multiplier, updated once per increment on acceptance (incremental Uzawa: lambda <- min(0, lambda + penalty * A * g)). The multiplier is constant within an increment (zero tangent contribution), drives the penetration toward zero over the increments at a fixed penalty, and sharpens the friction cone mu * N. Requires sliding=small. |
testfiles/edelweiss-only/NodeToDeformableSurfaceContact/test.inp*material, name=LinearElastic, id=linearelastic, provider=edelweiss
1.8e4, 0.3
*job, name=contactjob, domain=3d
*solver, solver=NIST, name=theSolver
*modelGenerator, generator=boxGen, name=lower
nX =3
nY =3
nZ =3
x0 =0
y0 =0
z0 =0
lX =2
lY =2
lZ =1
elProvider =edelweiss
elType =C3D8
*modelGenerator, generator=boxGen, name=upper
nX =3
nY =3
nZ =3
x0 =0.3
y0 =1.1
z0 =1.05
lX =0.6
lY =0.6
lZ =1
elProvider =edelweiss
elType =C3D8
*modelGenerator, generator=surfaceElementGenerator, name=gen1
surface = upper_back
name = slaveSurf
*modelGenerator, generator=surfaceElementGenerator, name=gen2
surface = lower_front
name = masterSurf
*section, name=section1, material=linearelastic, type=solid
lower_all
upper_all
*constraint, name=contact, type=nodeToDeformableSurfacePenalty
slaveSurface=slaveSurf_facets, masterSurface=masterSurf_facets, penalty=5e6, type=linear, searchDistance=2.0
*fieldOutput
>>perNode, elSet=lower_all, field=displacement, result=U, name=dispLower
>>perNode, elSet=upper_all, field=displacement, result=U, name=dispUpper
*output, type=ensight, name=esExport
>>perNode, fieldOutput=dispLower
>>perNode, fieldOutput=dispUpper
>>configuration, overwrite=yes
*step, solver=theSolver
maxInc=0.02, minInc=1e-6, maxNumInc=200, maxIter=25, stepLength=1
>>options, name=theSolver, extrapolation=off
>>dirichlet, name=fixLower, nSet=lower_back, field=displacement, 1=0.0, 2=0.0, 3=0.0
>>dirichlet, name=pushDown, nSet=upper_front, field=displacement, 3=-0.2
>>dirichlet, name=pinUpperXY, nSet=upper_bottomLeftBack, field=displacement, 1=0.0, 2=0.0
>>dirichlet, name=pinUpperRotZ, nSet=upper_bottomRightBack, field=displacement, 2=0.0
testfiles/edelweiss-only/NodeToDeformableSurfaceContactFrictionHexa20/test.inp*material, name=LinearElastic, id=linearelastic, provider=edelweiss
1.8e4, 0.3
*job, name=contactjob, domain=3d
*solver, solver=NIST, name=theSolver
*modelGenerator, generator=boxGen, name=lower
nX =2
nY =2
nZ =1
x0 =0
y0 =0
z0 =0
lX =2
lY =2
lZ =1
elProvider =edelweiss
elType =C3D20
*modelGenerator, generator=boxGen, name=upper
nX =2
nY =2
nZ =2
x0 =0.3
y0 =1.1
z0 =1.05
lX =0.6
lY =0.6
lZ =1
elProvider =edelweiss
elType =C3D20
*modelGenerator, generator=surfaceElementGenerator, name=gen1
surface = upper_back
name = slaveSurf
triangulation = midside
*modelGenerator, generator=surfaceElementGenerator, name=gen2
surface = lower_front
name = masterSurf
triangulation = midside
*section, name=section1, material=linearelastic, type=solid
lower_all
upper_all
** the quadratic force law is deliberate: its stiffness vanishes continuously at gap activation,
** which keeps the frictional tangent continuous for nodes lifting off/touching down at the
** tilting trailing/leading edges during the drag (the linear law's activation stiffness jump
** makes Newton limit-cycle there)
*constraint, name=contact, type=nodeToDeformableSurfacePenalty
slaveSurface=slaveSurf_facets, masterSurface=masterSurf_facets, penalty=3e7, type=quadratic,
searchDistance=2.0, sliding=small, mu=0.3, tangentPenalty=1e6
*fieldOutput
>>perNode, elSet=lower_all, field=displacement, result=U, name=dispLower
>>perNode, elSet=upper_all, field=displacement, result=U, name=dispUpper
>>perNode, name=RFDrag, nSet=lower_back, field=displacement, result=P, f(x)='sum(x[:,0])', saveHistory=True
>>perNode, name=RFNormal, nSet=lower_back, field=displacement, result=P, f(x)='sum(x[:,2])', saveHistory=True
*output, type=ensight, name=esExport
>>perNode, fieldOutput=dispLower
>>perNode, fieldOutput=dispUpper
>>configuration, overwrite=yes
*output, type=monitor, name=monitor
fieldOutput=RFDrag
fieldOutput=RFNormal
** step 1: press the small upper block onto the large lower block (closes the initial 0.05 gap,
** then compresses by 0.05)
*step, solver=theSolver
maxInc=0.34, minInc=1e-6, maxNumInc=200, maxIter=25, stepLength=1
>>options, name=theSolver, extrapolation=off
>>dirichlet, name=fixLower, nSet=lower_back, field=displacement, 1=0.0, 2=0.0, 3=0.0
>>dirichlet, name=holdUpper, nSet=upper_front, field=displacement, 1=0.0, 2=0.0, 3=-0.1
** step 2: drag the compressed upper block 0.4 across the lower block's surface -- the slave
** nodes slide over several master facet columns, exercising per-increment reassignment of the
** frozen small-sliding projection and the frictional stick-slip history
*step, solver=theSolver
maxInc=0.1, minInc=1e-6, maxNumInc=200, maxIter=25, stepLength=1
>>options, name=theSolver, extrapolation=off
>>dirichlet, name=fixLower, nSet=lower_back, field=displacement, 1=0.0, 2=0.0, 3=0.0
>>dirichlet, name=holdUpper, nSet=upper_front, field=displacement, 1=0.4, 2=0.0, 3=0.0
surfacetodeformablesurfacepenalty - Integrated surface-to-surface contact
The integrated (Gauss-point-to-segment) counterpart of the constraint above: contact is evaluated at quadrature points over the slave surface’s facets, and the resulting pressure is distributed with each side’s parent element face shape functions rather than with the flat facets’ own.
Use this instead of nodetodeformablesurfacepenalty on quadratic (hexa20/quad8) contact
surfaces. There, a node-based penalty scheme cannot transmit a correct pressure at all – the
consistent nodal load of a uniform pressure at a serendipity corner is tensile, which a
unilateral spring cannot exert, so the discrete solution opens every corner gap instead. The
mechanism, the measurements, and why no choice of nodal weights can repair it are documented in
Measured: the discrete solution lifts the corners off instead; the integrated formulation and its own limits are in
Integrated contact: distributing with the parent face’s basis. On matched linear faces the parent-face basis coincides with the facet
basis, so the two constraints agree to 13 digits and there is nothing to gain – but on
non-matching meshes the integrated formulation is a substantial improvement at any element
order, linear included, as the figure in Integrated contact: distributing with the parent face’s basis shows.
Scope: normal penalty contact under sliding=small, on both the implicit and the explicit solver
path. Coulomb friction and augmented Lagrange are not implemented here (the node-based
constraint has both), and sliding=finite is rejected rather than approximated.
Module edelweissfe.constraints.surfacetodeformablesurfacepenalty
Option |
Type / Default |
Description |
|---|---|---|
slaveSurface |
str (required) |
The element set of contact facet elements (Tria3ContactFacet/Line2ContactFacet) forming the slave surface; contact is integrated over these facets at quadrature points. The facets must come from the surface element generator, which stamps their parent faces. |
penalty |
float (required) |
The numerical penalty value, an interface stiffness modulus per unit slave surface area. |
type |
str, default=’linear’ |
The formulation type: ‘linear’ (linear force, constant stiffness with jump) or ‘quadratic’ (quadratic force, linear stiffness). |
nQuadraturePoints |
int, default=3 |
The number of quadrature points per slave facet: 1, 3 or 6 for a Tria3 facet (3D), 1, 2 or 3 for a Line2 facet (2D). The default of 3 integrates the parent face’s shape functions exactly over a facet, which is what makes the consistent nodal loads – including the negative corner loads of a serendipity face – come out exactly. |
searchDistance |
float, default=None |
An optional broadphase distance for the per-increment candidate-facet search. If not given, every contact point is always assigned its single closest master facet, without a distance gate. |
sliding |
str, default=’small’ |
The kinematic treatment of the contact geometry. Only ‘small’ (Abaqus-style small sliding: the closest-point projection – master facet, parametric location and normal – is frozen at every contact search from the last converged configuration, making the gap linear in the displacement DOFs) is implemented; ‘finite’ is rejected. |
masterSurface |
str (required) |
The element set of contact facet elements (Tria3ContactFacet/Line2ContactFacet) forming the master surface. |
testfiles/edelweiss-only/SurfaceToDeformableSurfaceContactPatchHexa20/test.inp** THE ACCEPTANCE TEST for integrated (Gauss-point-to-segment) contact on serendipity elements,
** i.e. for issue #112. Same model as NodeToDeformableSurfaceContactPatchHexa20 -- two blocks of
** identical 1x1 footprint and matching 2x2 hexa20 interface meshes under uniform compression,
** nu = 0 -- but with type=surfaceToDeformableSurfacePenalty. The exact solution is homogeneous
** uniaxial stress: contact pressure 360, total interface force 360, and a perfectly FLAT interface.
**
** The node-based constraint cannot represent this state at all. A quad8 face's consistent nodal
** loads under a uniform pressure are NEGATIVE at the corners (-pA/12, +pA/3 at the midsides), and a
** unilateral per-node spring cannot pull, so the discrete solution opens every corner gap instead:
** 0 of 9 corner nodes in contact, gaps of 2.4e-3, interface force 1.7% short, and a visibly
** undulating interface. Integrating the pressure and distributing it with the PARENT FACE's own
** quadratic shape functions removes the obstruction, because a corner's nodal force is then
** sum_q w_q p_q N_corner(xi_q), which is negative wherever N_corner is.
**
** Measured, at full load:
**
** force corners in contact max gap flatness
** C3D8 (linear), for reference 359.892 -- -1.2e-5 1.4e-17
** hexa20, nodeToDeformableSurface 353.787 0 of 9 +2.396e-3 1.4274e-3
** hexa20, surfaceToDeformableSurf 359.892 all 72 qp in contact -1.1996e-5 3.19e-16
**
** Every number is what theory demands, not merely "better":
** RFNormal 359.89203239028313 vs the C3D8 reference's 359.89203239028274 -- identical to 13
** digits, i.e. the same answer a linear interface gives, penalty compliance included.
** Uzspread 3.19e-16: the interface is flat to machine precision.
** gapMax -1.1996401079e-5 == -360/3e7 == -p/penalty at EVERY quadrature point: a uniform
** penetration, hence a uniform pressure. Nothing has lifted off.
** pNspread 5.55e-11, the patch test passing exactly, as its C3D8 sibling does.
** cornerForce (the minimum slave nodal normal force, negative = pulled) -29.99100269952982,
** against the analytic -4 * p * A / 12 = -29.99100269919026 for an interior corner
** shared by four faces of area A = 0.25. The corner nodes carry exactly the tensile
** load the face demands -- which is the entire point of the formulation.
**
** Note the pressures and gaps are per QUADRATURE POINT (24 slave facets x 3 points), so their
** fromExpression outputs are tied to the facet element set rather than to a node set; only
** cornerForce is nodal. See test_facetConsistent.inp of the node-based sibling for the baseline.
*material, name=LinearElastic, id=linearelastic, provider=edelweiss
1.8e4, 0.0
*job, name=patchjob, domain=3d
*solver, solver=NIST, name=theSolver
*modelGenerator, generator=boxGen, name=lower
nX =2
nY =2
nZ =1
x0 =0
y0 =0
z0 =0
lX =1
lY =1
lZ =1
elProvider =edelweiss
elType =C3D20
*modelGenerator, generator=boxGen, name=upper
nX =2
nY =2
nZ =1
x0 =0
y0 =0
z0 =1.02
lX =1
lY =1
lZ =1
elProvider =edelweiss
elType =C3D20
*modelGenerator, generator=surfaceElementGenerator, name=gen1
surface = upper_back
name = slaveSurf
triangulation = midside
*modelGenerator, generator=surfaceElementGenerator, name=gen2
surface = lower_front
name = masterSurf
triangulation = midside
*section, name=section1, material=linearelastic, type=solid
lower_all
upper_all
*constraint, name=contact, type=surfaceToDeformableSurfacePenalty
slaveSurface=slaveSurf_facets, masterSurface=masterSurf_facets, penalty=3e7, type=linear,
searchDistance=1.0, sliding=small, nQuadraturePoints=3
*fieldOutput
>>perNode, elSet=lower_all, field=displacement, result=U, name=dispLower
>>perNode, elSet=upper_all, field=displacement, result=U, name=dispUpper
>>perNode, name=RFNormal, nSet=lower_back, field=displacement, result=P, f(x)='sum(x[:,2])', saveHistory=True
>>perNode, name=Uzspread, nSet=slaveSurf_nodes, field=displacement, result=U, f(x)='np.max(x[:,2])-np.min(x[:,2])', saveHistory=True
>>fromExpression, name=gapMax, elSet=slaveSurf_facets, expression='model.constraints["contact"].getGaps()', f(x)='np.max(x)', saveHistory=True
>>fromExpression, name=cornerForce, nSet=slaveSurf_nodes, expression='model.constraints["contact"].getSlaveNodalNormalForces()', f(x)='np.min(x)', saveHistory=True
>>fromExpression, name=pNmean, elSet=slaveSurf_facets, expression='model.constraints["contact"].getNormalPressures()', f(x)='np.mean(x[x>0]) if np.any(x>0) else np.float64(0.0)', saveHistory=True
>>fromExpression, name=pNspread, elSet=slaveSurf_facets, expression='model.constraints["contact"].getNormalPressures()', f(x)='(np.max(x[x>0])-np.min(x[x>0]))/np.mean(x[x>0]) if np.any(x>0) else np.float64(0.0)', saveHistory=True
*output, type=monitor, name=monitor
fieldOutput=RFNormal
fieldOutput=Uzspread
fieldOutput=gapMax
fieldOutput=cornerForce
fieldOutput=pNmean
fieldOutput=pNspread
*step, solver=theSolver
maxInc=0.25, minInc=1e-6, maxNumInc=100, maxIter=25, stepLength=1
>>options, name=theSolver, extrapolation=off
>>dirichlet, name=fixLower, nSet=lower_back, field=displacement, 1=0.0, 2=0.0, 3=0.0
>>dirichlet, name=pushUpper, nSet=upper_front, field=displacement, 1=0.0, 2=0.0, 3=-0.06
testfiles/edelweiss-only/NEDSurfaceContact/test.inp** Integrated (Gauss-point-to-segment) contact on the EXPLICIT dynamic path: the same model and the
** same decisive check as NEDContact, with type=surfaceToDeformableSurfacePenalty.
**
** As there, the check is not that this runs but the `dispLowerMax` history: nothing else in this
** model loads the lower block, so its displacement is identically zero unless the contact
** constraint's forces actually reach the equations. `cornerForce` is the minimum slave nodal normal
** force; it is POSITIVE here because these are linear (C3D8) faces, whose consistent nodal loads
** are all non-negative -- contrast SurfaceToDeformableSurfaceContactPatchHexa20, where the same
** quantity is negative by design.
**
** DETERMINISM, and why contact-update-frequency=0: the closest-point projection is established
** once before the increment loop and never re-searched. That is a requirement, not a convenience.
** With a periodic re-search this model cannot serve as a bitwise reference at all -- perturbing the
** penalty by 2e-16 relative (5e4 -> 50000.00000000001) and running to completion moves the final
** displacements by 5.3e-4 at contact-update-frequency=10, against the runner's fixed 1e-6 absolute
** tolerance. Pinned, the same perturbation ends at 3.8e-16 and the case reproduces across thread
** count and machine. Pinning costs about 1% in the monitored results (dispLowerMax 0.03877 against
** 0.03842), which is the small-sliding approximation error this option exists to trade away, not a
** loss of coverage: the constraint's force path, its explicit hook and its active set are exercised
** exactly as before.
**
** The switch is master-facet REASSIGNMENT, and it is worth stating precisely because the obvious
** suspects are both innocent. It is not the `g < 0` activation, and it is not the batched
** evaluation's summation order (that path agrees with the per-point reference loop to 2e-15
** relative, and exactly for the linear law). Over this run 166 contact points cross from one master
** element's parent face onto the next one's; at a crossing the frozen facet normal changes by up to
** 5.2e-2 -- the kink between two adjacent, now non-coplanar, deformed element faces -- which jumps
** that point's gap by up to 4.8e-4, some 1.3% of the gap itself. Which side of the seam a
** borderline point lands on is decided at the 1e-16 level, so the trajectory is a step function of
** rounding.
**
** Re-projection by itself is NOT the problem, and pinning is not a way of hiding one. The
** node-based constraint under sliding=small injects a gap jump at every update once contact
** engages (up to 4.8e-6 on the equivalent model) and still reproduces to 6.9e-16, because its 16
** slave nodes never cross a master element boundary, while 54 quadrature points spread over the
** slave surface do. The discrete facet switch is the chaotic ingredient; the continuous re-freeze
** is not. Note also that NodeToDeformableSurfaceContact-style decks do not set `sliding` at all and
** so run the node-based default, `finite`, which has no frozen projection to re-freeze -- they are
** not a control for any of this unless sliding=small is set explicitly.
**
** What this deck therefore does not cover is updateConnectivity's reassignment path, which is
** covered exactly, and deterministically, in test_surfacetodeformablesurfacepenalty.py:
** re-projection within one parent face, and across two parent faces of a flat master, are both
** asserted to leave the gap bit-identical. Do not reintroduce a periodic update here to win that
** coverage back -- it buys a flaky reference for a mechanism a unit test already pins exactly.
**
** Deliberately C3D8 rather than C3D20, although hexa20 is the reason the constraint exists: the
** pure-Python C3D20 element cannot currently run explicit dynamics at all. Row-sum lumping of its
** consistent mass matrix yields NEGATIVE corner masses (verified pre-existing: unmodified
** NEDContact with only elType changed to C3D20 fails identically, with the node-based contact),
** which is the very same serendipity property that motivates this constraint, appearing in the
** mass matrix instead of the contact loads -- and it wants the same kind of fix there
** (Duczek & Gravenkamp's 1/3-2/3 lumping weights). Marmot's own quadratic elements lump their mass
** themselves and are unaffected, so this does not block hexa20 explicit models that use them.
*material, name=LinearElastic, id=linearelastic, provider=edelweiss
1.8e4, 0.3, 1.0
*job, name=nedcontactjob, domain=3d
*solver, solver=NEDParallel, name=theSolver
courant-number=0.8
output-frequency=50
contact-update-frequency=0
*modelGenerator, generator=boxGen, name=lower
nX =3
nY =3
nZ =3
x0 =0
y0 =0
z0 =0
lX =2
lY =2
lZ =1
elProvider =edelweiss
elType =C3D8
*modelGenerator, generator=boxGen, name=upper
nX =3
nY =3
nZ =3
x0 =0.3
y0 =1.1
z0 =1.05
lX =0.6
lY =0.6
lZ =1
elProvider =edelweiss
elType =C3D8
*modelGenerator, generator=surfaceElementGenerator, name=gen1
surface = upper_back
name = slaveSurf
*modelGenerator, generator=surfaceElementGenerator, name=gen2
surface = lower_front
name = masterSurf
*section, name=section1, material=linearelastic, type=solid
lower_all
upper_all
*constraint, name=contact, type=surfaceToDeformableSurfacePenalty
slaveSurface=slaveSurf_facets, masterSurface=masterSurf_facets, penalty=5e4, type=linear, searchDistance=2.0
*fieldOutput
>>perNode, elSet=lower_all, field=displacement, result=U, name=dispLower
>>perNode, elSet=upper_all, field=displacement, result=U, name=dispUpper
>>perNode, name=dispLowerMax, elSet=lower_all, field=displacement, result=U, f(x)='np.max(np.abs(x))', saveHistory=True
>>perNode, name=RFLower, nSet=lower_back, field=displacement, result=P, f(x)='np.sum(x[:,2])', saveHistory=True
>>fromExpression, name=cornerForce, nSet=slaveSurf_nodes, expression='model.constraints["contact"].getSlaveNodalNormalForces()', f(x)='np.min(x)', saveHistory=True
*output, type=monitor, name=monitor
fieldOutput=dispLowerMax
fieldOutput=RFLower
fieldOutput=cornerForce
*step, type=adaptiveForExplicitSimulations, solver=theSolver
maxInc=1, minInc=1e-12, maxNumInc=2000, maxIter=25, stepLength=0.5
>>dirichlet, name=fixLower, nSet=lower_back, field=displacement, 1=0.0, 2=0.0, 3=0.0
>>dirichlet, name=pushDown, nSet=upper_front, field=displacement, 3=-0.2
>>dirichlet, name=pinUpperXY, nSet=upper_bottomLeftBack, field=displacement, 1=0.0, 2=0.0
>>dirichlet, name=pinUpperRotZ, nSet=upper_bottomRightBack, field=displacement, 2=0.0
surfacetodiscreterigidbodypenalty - Integrated contact against a discrete rigid body
The integrated counterpart of nodetodiscreterigidbodypenalty, and the rigid-body counterpart
of surfacetodeformablesurfacepenalty: contact is evaluated at quadrature points over a slave
facet surface and distributed with the slave parent-face shape functions, against the triangulated,
closed surface of a discrete rigid body. The theory is documented in
Integrated contact against a discrete rigid body. Use it for any contact against rigid supports or indenters,
in particular with quadratic slave faces.
Scope: normal penalty contact under sliding=small, 3D, implicit and explicit solver paths, live
AMR of the slave side. The slave surface must consist of a single element type.
Module edelweissfe.constraints.surfacetodiscreterigidbodypenalty
Option |
Type / Default |
Description |
|---|---|---|
slaveSurface |
str (required) |
The element set of contact facet elements (Tria3ContactFacet/Line2ContactFacet) forming the slave surface; contact is integrated over these facets at quadrature points. The facets must come from the surface element generator, which stamps their parent faces. |
penalty |
float (required) |
The numerical penalty value, an interface stiffness modulus per unit slave surface area. |
type |
str, default=’linear’ |
The formulation type: ‘linear’ (linear force, constant stiffness with jump) or ‘quadratic’ (quadratic force, linear stiffness). |
nQuadraturePoints |
int, default=3 |
The number of quadrature points per slave facet: 1, 3 or 6 for a Tria3 facet (3D), 1, 2 or 3 for a Line2 facet (2D). The default of 3 integrates the parent face’s shape functions exactly over a facet, which is what makes the consistent nodal loads – including the negative corner loads of a serendipity face – come out exactly. |
searchDistance |
float, default=None |
An optional broadphase distance for the per-increment candidate-facet search. If not given, every contact point is always assigned its single closest master facet, without a distance gate. |
sliding |
str, default=’small’ |
The kinematic treatment of the contact geometry. Only ‘small’ (Abaqus-style small sliding: the closest-point projection – master facet, parametric location and normal – is frozen at every contact search from the last converged configuration, making the gap linear in the displacement DOFs) is implemented; ‘finite’ is rejected. |
rigidBody |
str (required) |
The name of the discrete rigid body (as registered in model.rigidBodies). Its surface must be closed and consist of triangles. |
testfiles/edelweiss-only/SurfaceToDiscreteRigidBodyContactPatchHexa20/test.inp** Integrated (Gauss-point-to-segment) contact of a serendipity (hexa20) block against a DISCRETE
** RIGID BODY. The block (1x1 footprint, 2x2x1 hexa20, nu = 0) is pressed onto a fixed rigid box
** (an STL, larger than the footprint) whose top face is at z = 1. The exact solution is homogeneous
** uniaxial stress with a perfectly FLAT interface.
**
** The node-based nodeToDiscreteRigidBodyPenalty cannot represent this state: a quad8 face's
** consistent nodal loads under a uniform pressure are NEGATIVE at the corners (-pA/12, +pA/3 at the
** midsides), which a unilateral per-node spring cannot exert, so the corners lift off. Here the
** pressure is integrated at quadrature points and distributed with the PARENT FACE's quadratic
** shape functions, so the corner nodes carry the tensile load the face demands.
**
** Checked quantities: cornerForce (the minimum slave nodal normal force, negative = pulled) equals
** -4 p A / 12 for an interior corner shared by four faces of area A = 0.25; the interface is flat
** (Uzspread); every quadrature point is closed with the same gap -p/penalty (gapMax); the pressure
** is uniform (pNspread); and the total contact force (contactForce, negative in compression)
** balances the reaction at the loaded face (RFblock).
*material, name=LinearElastic, id=linearelastic, provider=edelweiss
1.8e4, 0.0
*job, name=patchjob, domain=3d
*solver, solver=NIST, name=theSolver
*modelGenerator, generator=boxGen, name=block
nX =2
nY =2
nZ =1
x0 =0
y0 =0
z0 =1.02
lX =1
lY =1
lZ =1
elProvider =edelweiss
elType =C3D20
*modelGenerator, generator=surfaceElementGenerator, name=gen1
surface = block_back
name = slaveSurf
triangulation = midside
*modelGenerator, generator=discreteRigidBodyGenerator, name=support, executeAfterManualGeneration=True
filename=support_box.stl
rpCoordinate='0.5, 0.5, 0.5'
*section, name=section1, material=linearelastic, type=solid
block_all
*constraint, name=contact, type=surfaceToDiscreteRigidBodyPenalty
slaveSurface=slaveSurf_facets, rigidBody=support, penalty=3e7, type=linear,
searchDistance=1.0, nQuadraturePoints=3
*fieldOutput
>>perNode, elSet=block_all, field=displacement, result=U, name=displacement
>>perNode, name=RFblock, nSet=block_front, field=displacement, result=P, f(x)='sum(x[:,2])', saveHistory=True
>>fromExpression, name=contactForce, nSet=support_rp, expression='np.array([model.constraints["contact"].totalNormalForce])', f(x)='x[0]', saveHistory=True
>>perNode, name=Uzspread, nSet=slaveSurf_nodes, field=displacement, result=U, f(x)='np.max(x[:,2])-np.min(x[:,2])', saveHistory=True
>>fromExpression, name=gapMax, elSet=slaveSurf_facets, expression='model.constraints["contact"].getGaps()', f(x)='np.max(x)', saveHistory=True
>>fromExpression, name=cornerForce, nSet=slaveSurf_nodes, expression='model.constraints["contact"].getSlaveNodalNormalForces()', f(x)='np.min(x)', saveHistory=True
>>fromExpression, name=pNmean, elSet=slaveSurf_facets, expression='model.constraints["contact"].getNormalPressures()', f(x)='np.mean(x[x>0]) if np.any(x>0) else np.float64(0.0)', saveHistory=True
>>fromExpression, name=pNspread, elSet=slaveSurf_facets, expression='model.constraints["contact"].getNormalPressures()', f(x)='(np.max(x[x>0])-np.min(x[x>0]))/np.mean(x[x>0]) if np.any(x>0) else np.float64(0.0)', saveHistory=True
*output, type=monitor, name=monitor
fieldOutput=RFblock
fieldOutput=contactForce
fieldOutput=Uzspread
fieldOutput=gapMax
fieldOutput=cornerForce
fieldOutput=pNmean
fieldOutput=pNspread
*step, solver=theSolver
maxInc=0.25, minInc=1e-6, maxNumInc=100, maxIter=25, stepLength=1
>>options, name=theSolver, extrapolation=off
>>dirichlet, name=fixSupportDisplacement, nSet=support_rp, field=displacement, 1=0.0, 2=0.0, 3=0.0
>>dirichlet, name=fixSupportRotation, nSet=support_rp, field=rotation, 1=0.0, 2=0.0, 3=0.0
>>dirichlet, name=pushBlock, nSet=block_front, field=displacement, 1=0.0, 2=0.0, 3=-0.06
testfiles/edelweiss-only/NEDSurfaceToDiscreteRigidBodyContact/test.inp** Integrated contact against a DISCRETE RIGID BODY on the EXPLICIT dynamic path: a block is pushed
** onto a fixed rigid box (an STL, top face at z = 1), with the same material, loading and solver
** settings as NEDSurfaceContact, whose deformable lower block is here replaced by the rigid box.
**
** The decisive checks: penetrationMax, the deepest penetration of the block into the rigid top
** face, stays at the penalty compliance (about 0.04 = pressure / penalty at the end) instead of
** growing with the imposed push of 0.2; and contactForce, the total normal force of the constraint,
** follows the reaction at the pushed face (RFUpper) up to the inertia of the block.
**
** contact-update-frequency=10 exercises the periodic contact search. Unlike NEDSurfaceContact,
** this stays deterministic: the rigid top face is flat, so a contact point moving from one of its
** triangles onto the next sees the same plane, and its gap does not jump.
**
** The explicit solver needs inertia on every degree of freedom of a field, so the rigid body is
** given a (here irrelevant, since its reference point is fixed) mass and rotational inertia.
**
** C3D8 rather than C3D20, because the pure-Python C3D20 element cannot run explicit dynamics (its
** row-sum lumped mass has negative corner entries), see NEDSurfaceContact.
*material, name=LinearElastic, id=linearelastic, provider=edelweiss
1.8e4, 0.3, 1.0
*job, name=nedrigidcontactjob, domain=3d
*solver, solver=NEDParallel, name=theSolver
courant-number=0.8
output-frequency=50
contact-update-frequency=10
*modelGenerator, generator=boxGen, name=upper
nX =3
nY =3
nZ =3
x0 =0.3
y0 =1.1
z0 =1.05
lX =0.6
lY =0.6
lZ =1
elProvider =edelweiss
elType =C3D8
*modelGenerator, generator=surfaceElementGenerator, name=gen1
surface = upper_back
name = slaveSurf
*modelGenerator, generator=discreteRigidBodyGenerator, name=support, executeAfterManualGeneration=True
filename=support_box.stl
rpCoordinate='1.0, 1.0, 0.5'
mass=1.0
inertia="1.0, 1.0, 1.0"
*section, name=section1, material=linearelastic, type=solid
upper_all
*constraint, name=contact, type=surfaceToDiscreteRigidBodyPenalty
slaveSurface=slaveSurf_facets, rigidBody=support, penalty=5e4, type=linear, searchDistance=2.0
*fieldOutput
>>perNode, elSet=upper_all, field=displacement, result=U, name=dispUpper
>>perNode, name=penetrationMax, nSet=slaveSurf_nodes, field=displacement, result=U, f(x)='np.maximum(0.0, -np.min(x[:,2]) - 0.05)', saveHistory=True
>>perNode, name=RFUpper, nSet=upper_front, field=displacement, result=P, f(x)='np.sum(x[:,2])', saveHistory=True
>>fromExpression, name=contactForce, nSet=support_rp, expression='np.array([model.constraints["contact"].totalNormalForce])', f(x)='x[0]', saveHistory=True
*output, type=monitor, name=monitor
fieldOutput=penetrationMax
fieldOutput=RFUpper
fieldOutput=contactForce
*step, type=adaptiveForExplicitSimulations, solver=theSolver
maxInc=1, minInc=1e-12, maxNumInc=2000, maxIter=25, stepLength=0.5
>>dirichlet, name=fixSupportDisplacement, nSet=support_rp, field=displacement, 1=0.0, 2=0.0, 3=0.0
>>dirichlet, name=fixSupportRotation, nSet=support_rp, field=rotation, 1=0.0, 2=0.0, 3=0.0
>>dirichlet, name=pushDown, nSet=upper_front, field=displacement, 3=-0.2
>>dirichlet, name=pinUpperXY, nSet=upper_bottomLeftBack, field=displacement, 1=0.0, 2=0.0
>>dirichlet, name=pinUpperRotZ, nSet=upper_bottomRightBack, field=displacement, 2=0.0
tie - Surface-to-surface tie (DOF elimination)
An Abaqus-style tie constraint bonding a slave surface rigidly to a deformable master surface.
Unlike the penalty- and Lagrange-multiplier-based constraints above, the tie is enforced by
master-slave DOF elimination (multi-point constraint condensation): each slave node’s
displacement is constrained to the frozen facet-linear interpolation of its closest master facet,
\(u_s = \sum_a N_a \, u_{m_a}\), and the slave DOFs are condensed out of the equation system
– exact enforcement, no penalty parameter, no additional DOFs, zero added stiffness. The same
mechanism serves the implicit solvers (system matrix transformation
\(\tilde{K} = T^T K \, T + I_s\), the identity \(I_s\) filling the structurally zero
slave rows and columns, followed by the recovery \(\delta u_s \mathrel{+}= \sum_a N_a \,
\delta u_{m_a}\) after each linear solve – symmetric whenever \(K\) is, and free of any scale of
its own, so it equilibrates cleanly also for a consistent mass matrix; see
MultiPointConstraintTransformation) and explicit
dynamics (mass-conserving row-sum folding of slave masses and forces onto the masters, direct
kinematic slaving – the critical time step is untouched).
The implicit solver’s useThreadedMPCCondensation option (Solvers, NIST) computes
\(T^T K T + I_s\) in place, with an OpenMP-threaded kernel
(MPCCondensation), instead of with SciPy’s
single-threaded CSR routines, which build a second matrix: the system matrix is held once, not twice.
For this, the CSR pattern of \(K\) gains, once per equation system, the entries where a master
collects the couplings of its slaves
(condensedPatternFillIn());
the slave columns of \(K\) become stored zeros in turn. Only the rows a constraint touches are
computed, each by one thread from the rows of \(K\) it combines, in a fixed order of summation, so
the result does not depend on the number of threads. Explicit zeros are pruned from a copy according to
pruneCondensedMatrixZeros. It is off by default.
On matching interface meshes, tying is exactly equivalent to merging the interface nodes: the
patch test passes to machine precision, and both the implicit and the explicit (central
difference) solutions reproduce the monolithic mesh identically. On non-matching meshes the
constraint is still enforced exactly, but the classical node-to-surface limitation applies: the
slave surface’s consistent nodal force pattern, redistributed by the facet-linear weights, does
not match the master’s own consistent pattern, so a uniform stress state exhibits a bounded,
refinement-convergent interface perturbation (quantified in
testfiles/edelweiss-only/TieHexa20Patch/test_mismatched.inp) – the same obstruction that
motivates mortar methods.
For quadratic (hexa20/quad8) faces, generate both facet surfaces with triangulation=midside:
the corner triangulation excludes the midside nodes from the facet node lists entirely, leaving
slave midside nodes untied.
Module edelweissfe.constraints.tie
Option |
Type / Default |
Description |
|---|---|---|
slaveSurface |
str (required) |
The element set of contact facet elements (Tria3ContactFacet/Line2ContactFacet) forming the slave surface; its nodes are tied. For quadratic (hexa20/quad8) faces, generate the facets with triangulation=midside on BOTH surfaces – the corner triangulation excludes the midside nodes from the facet node list entirely, leaving them untied. |
masterSurface |
str (required) |
The element set of contact facet elements (Tria3ContactFacet/Line2ContactFacet) forming the master surface. |
positionTolerance |
float, default=None |
Whether a slave node is tied at all: nodes whose reference-configuration closest-point distance to the master surface exceeds this tolerance are left untied (recorded in the constraint’s untiedSlaveNodes), matching Abaqus’ *TIE default behavior of silently dropping out-of-range slave nodes. Applies independently of adjust – whether a tied node’s coordinates additionally get snapped onto the master is a separate decision, see adjust/adjustTolerance. If not given (the default), a tolerance is computed as positionToleranceFactor times the master surface’s characteristic facet size – see positionToleranceFactor. Set this explicitly to an absolute distance to override that. |
positionToleranceFactor |
float, default=0.25 |
Used only when positionTolerance is not given: the default tolerance is this fraction of the master surface’s characteristic (mean, over all its facets) edge length, computed once and used for every slave node. 0.25 comfortably exceeds the sub-percent gaps expected between two compatible discretizations of the same surface (mismatched density, curvature/interpolation error), while remaining well below the facet-size-or-larger gaps that indicate the surfaces don’t actually correspond (e.g. a slave surface extending beyond the master surface’s actual extent – a partial-bond-length or otherwise partially-overlapping pair of surfaces). |
adjust |
bool, default=True |
Whether a TIED node’s coordinates additionally get snapped onto its projected master point at construction (Abaqus-like default) – a separate decision from whether the node is tied at all (see positionTolerance). If False, no tied node is ever snapped: any initial geometric gap is preserved rigidly (the displacements are tied, not the positions), regardless of size. If True, a tied node is snapped only if its distance is also within adjustTolerance (default: unconditionally, i.e. every tied node is snapped, matching plain Abaqus ADJUST=YES) – see adjustTolerance to snap away only small, effectively-numerical gaps while still tying (without snapping) across larger, deliberate ones. Note that adjusting modifies the nodal coordinates before the element geometry is initialized; avoid adjusting nodes that also belong to an already-generated contact surface of another constraint. |
adjustTolerance |
float, default=None |
Used only when adjust=True: a tied node’s coordinates are snapped onto the master only if its closest-point distance is also within this tolerance; beyond it, the node stays tied (kinematically) but its position is left as found, preserving the gap. Independent of positionTolerance – a node can be tied across a fairly generous distance while only genuinely small (e.g. sub-percent, mesh-discretization-scale) gaps within that get snapped away. If not given (the default), any tied node is snapped, matching plain Abaqus ADJUST=YES and this constraint’s behavior before this option existed. |
testfiles/edelweiss-only/TieHexa20Patch/test.inp** tie patch test: two hexa20 blocks with matching interface meshes (2x2 vs 2x2), bonded by the
** DOF-elimination tie constraint, under uniform compression with nu = 0. The upper block, 0.98
** tall, starts with a 0.02 gap that the tie's adjust (default) closes at construction by moving
** its bottom nodes -- and its elements with them -- onto the lower block, so the snapped column is
** exactly 2 tall and the exact solution is homogeneous uniaxial stress: sigma = E * 0.06 / 2 = 540
** (before 2026-10-02 the elements kept their pre-snap geometry, and a 1-tall upper block gave the
** same number for the wrong reason), and the interface displacement is
** perfectly planar. On matching meshes every slave node projects exactly onto a master node
** (weights snap to a facet vertex), so the patch test passes to machine precision: RFBottom =
** 540 exactly, and the spread of the slave nodes' z-displacement (UzSpread) is ~1e-17. See
** test_mismatched.inp for the non-matching variant quantifying the intrinsic node-to-surface tie
** limitation. triangulation=midside on both surfaces is mandatory for hexa20: it puts the
** midside nodes into the facet node lists (slave side: they must be tied; master side: they
** carry interpolation weights).
*material, name=LinearElastic, id=linearelastic, provider=edelweiss
1.8e4, 0.0
*job, name=tiejob, domain=3d
*solver, solver=NIST, name=theSolver
*modelGenerator, generator=boxGen, name=lower
nX =2
nY =2
nZ =1
x0 =0
y0 =0
z0 =0
lX =1
lY =1
lZ =1
elProvider =edelweiss
elType =C3D20
*modelGenerator, generator=boxGen, name=upper
nX =2
nY =2
nZ =1
x0 =0
y0 =0
z0 =1.02
lX =1
lY =1
lZ =0.98
elProvider =edelweiss
elType =C3D20
*modelGenerator, generator=surfaceElementGenerator, name=gen1
surface = upper_back
name = slaveSurf
triangulation = midside
*modelGenerator, generator=surfaceElementGenerator, name=gen2
surface = lower_front
name = masterSurf
triangulation = midside
*section, name=section1, material=linearelastic, type=solid
lower_all
upper_all
*constraint, name=tie1, type=tie
slaveSurface=slaveSurf_facets, masterSurface=masterSurf_facets
*fieldOutput
>>perNode, elSet=lower_all, field=displacement, result=U, name=dispLower
>>perNode, elSet=upper_all, field=displacement, result=U, name=dispUpper
>>perNode, name=RFBottom, nSet=lower_back, field=displacement, result=P, f(x)='sum(x[:,2])', saveHistory=True
>>perNode, name=UzSpread, nSet=slaveSurf_nodes, field=displacement, result=U, f(x)='np.max(x[:,2])-np.min(x[:,2])', saveHistory=True
*output, type=monitor, name=monitor
fieldOutput=RFBottom
fieldOutput=UzSpread
*step, solver=theSolver
maxInc=0.5, minInc=1e-6, maxNumInc=100, maxIter=25, stepLength=1
>>options, name=theSolver, extrapolation=off
>>dirichlet, name=fixBottom, nSet=lower_back, field=displacement, 1=0.0, 2=0.0, 3=0.0
>>dirichlet, name=pushTop, nSet=upper_front, field=displacement, 3=-0.06
testfiles/edelweiss-only/TieNED/test.inp** explicit-dynamics (NED, central difference) tie test: two CPE4 patches with MISMATCHED
** interface meshes (2 vs 3 elements across), bonded by the DOF-elimination tie, top edge pushed
** down by a prescribed-velocity Dirichlet ramp. The tie in explicit dynamics adds no stiffness
** (critical time step untouched): slave masses and forces are folded onto the masters
** (row-sum lumping of T^T M T, mass-conserving), and the slave kinematics are assigned directly
** from the masters each increment -- a compression wave crosses the tied interface. Regression
** via U.ref; the journal's internal/kinetic energy split provides the sanity monitor. See
** test_merged.inp for the manually-run monolithic reference of the same geometry.
*material, name=LinearElastic, id=linearelastic, provider=edelweiss
100, 0.0, 1.0
*job, name=tiejobned, domain=2d
*solver, solver=NED, name=theSolver
output-frequency=10
*modelGenerator, generator=planeRectQuad, name=lower
nX =2
nY =2
x0 =0
y0 =0
l =1
h =1
elType =CPE4
elProvider =edelweiss
*modelGenerator, generator=planeRectQuad, name=upper
nX =3
nY =2
x0 =0
y0 =1.02
l =1
h =1
elType =CPE4
elProvider =edelweiss
*modelGenerator, generator=surfaceElementGenerator, name=gen1
surface = upper_bottom
name = slaveSurf
*modelGenerator, generator=surfaceElementGenerator, name=gen2
surface = lower_top
name = masterSurf
*section, name=section1, thickness=1.0, material=linearelastic, type=plane
lower_all
upper_all
*constraint, name=tie1, type=tie
slaveSurface=slaveSurf_facets, masterSurface=masterSurf_facets
*fieldOutput
>>perNode, elSet=lower_all, field=displacement, result=U, name=dispLower
>>perNode, elSet=upper_all, field=displacement, result=U, name=dispUpper
>>perNode, name=RFBottom, nSet=lower_bottom, field=displacement, result=P, f(x)='sum(x[:,1])', saveHistory=True
*output, type=monitor, name=monitor
fieldOutput=RFBottom
*step, type=adaptiveForExplicitSimulations, solver=theSolver
maxInc=1e0, minInc=1e-2, maxNumInc=10000, maxIter=25, stepLength=1
>>dirichlet, name=fixBottom, nSet=lower_bottom, field=displacement, 1=0, 2=0
>>dirichlet, name=pushTop, nSet=upper_top, field=displacement, 2=-0.05
Contested slave degrees of freedom
A degree of freedom can be condensed out only once, so when several multi-point constraints ask for
the same slave DOF, exactly one may keep it. This arises routinely under \(h\)-adaptivity: a
hanging node created on a refined tie interface is a hangingnode slave, and the tie would like to
constrain it too.
The decision is made centrally, where every constraint’s records are collected, and not inside any constraint: the first constraint in model order keeps the DOF, and later claims on it are dropped and reported. Resolving it centrally is what makes the result independent of the order in which constraints are refreshed – a constraint that inspected its peers mid-refresh would read some of them in their pre-refinement state.
Model order therefore carries meaning. For the hanging-node case the precedence is not
arbitrary, and it is expressed by hAdaptivity registering its hanging-node constraint at the
front of the model’s multi-point constraints:
If the hanging-node constraint keeps the node, the refined and unrefined meshes stay exactly conforming, and the node is still tied – indirectly but almost exactly, because its coarse parent nodes are themselves tie slaves and it is their interpolation.
If the tie keeps it instead, the node is bonded exactly, but nothing holds it on its coarse parent edge any more, and the refined and unrefined meshes come apart there.
One constraint has a stand-in for its condition and the other does not, which is what settles it. Measured on the coarse anchor pry-out, the two precedences differ by 5.3e-02 relative displacement and eventually by the refined mesh itself.
Note
This used to be resolved inside the tie, which built its records against its peers’
claimedSlaveNodes() while the refresh sweep was still running – so peers that had not yet
refreshed answered with their pre-refinement claims. A node that had been a hanging node before
the latest refinement, and was not one after it, was reported as still claimed, and the tie
dropped it, leaving it constrained by nothing at all.
Multi-point constraints and Dirichlet boundary conditions
A degree of freedom eliminated by a multi-point constraint – a tie slave, a hangingnode
slave – cannot also be prescribed by a Dirichlet boundary condition: the constraint already
determines it from its masters, and a prescribed value would be a second, independent equation for
the same unknown. Where a boundary condition’s node set overlaps a constraint’s slave nodes, the
solvers resolve the overlap automatically: the conflicting constraint equation is dropped,
for that component only, so that the boundary condition takes precedence. This is how Abaqus
resolves the same overlap, and it needs no configuration.
The alternative – rejecting such a model and requiring the user to subtract the constraint’s slave nodes from the boundary condition’s node set – is worth understanding, because input files written for solvers that demand it often carry such a subtraction. It is a snapshot: computed against one mesh, one tie tolerance and one refinement state, and silently no longer correct when any of them changes, with no diagnostic, because what remains is a plain list of node labels that no longer records why those labels were chosen. It also interferes with \(h\)-adaptivity, which extends a node set to newly created nodes only where the whole parent face or edge already lies within the set: every face touching a subtracted node stops propagating the boundary condition to its children, so the boundary condition thins out with each refinement. Such subtractions should be removed rather than maintained.
Dropped equations are reported in two classes, and the distinction matters:
- redundant
Every master carrying non-negligible weight is itself prescribed, and the weighted sum reproduces the value prescribed on the slave. The constraint equation and the boundary condition say the same thing, so dropping either changes nothing. Reported as a count, at journal level 1. This is the expected case for a symmetry plane cutting through a tied interface: the slave and its masters all lie on the plane and all carry the same prescribed value.
- not implied
They do not say the same thing – some master is unprescribed, or the weighted sum disagrees. The boundary condition still wins, but the model has genuinely changed, so this is reported as a warning naming how much of the constraint weight rests on unprescribed masters. It is a signal worth chasing rather than silencing: the usual cause is a boundary condition whose node set is incomplete, so that masters which ought to carry the condition do not.
Masters are resolved transitively before classification – a master may itself be a slave of another constraint – and weights are compared against a scaled tolerance, never against exact zero, since clamped closest-point projections routinely leave weights at round-off magnitude.
hangingnode - Hanging-node coupling for adaptive mesh refinement
The kinematic multi-point constraint that couples a non-conforming 2:1 refined interface produced by hanging-node \(h\)-adaptivity (see Model modifiers). Each hanging (slave) node is tied to the coarse serendipity trace it lies on – an 8-node QUAD8 face or a 3-node quadratic edge – by \(u_s = \sum_a N_a(\xi_s)\, u_{m_a}\), enforced by master-slave DOF elimination (no Lagrange multipliers, no extra DOFs, no saddle point). Because the QUAD8 face-trace and quadratic edge spaces are nested under octree refinement, the coupling is exact. All fields active on the node (displacement, nonlocal damage, …) are constrained with the same field-independent weights. Multi-level chains are pre-flattened to independent masters, as required by the DOF-elimination transformation.
The records (masters + weights per slave) are either loaded from a file (recordsFile, one line
per slave <slaveLabel> <masterLabel> <weight> ...) for a statically pre-refined mesh, or set in
memory by the hAdaptivity model modifier after each dynamic refinement.
Module edelweissfe.constraints.hangingnode
Hanging-node constraint for adaptive-mesh-refinement, enforced by master-slave DOF elimination (multi-point constraint), Abaqus-style – no Lagrange multipliers, no extra DOFs, no saddle point.
Each hanging (slave) node on a 2:1-refined HEX20 face is tied to the coarse serendipity trace it hangs on:
for EVERY field active on the node (displacement, nonlocal damage, …) and every component. Because
the QUAD8 face-trace / quadratic-edge spaces are nested under octree refinement, these constraints
are exact. The records are precomputed (masters + weights, with multi-level chains already flattened
to independent masters) by the adaptivity machinery – see
edelweissfe.adaptivity.refinement.AdaptiveMesh.hanging_mpc_records() – and written to a plain
records file, one line per slave node:
<slaveLabel> <masterLabel> <weight> <masterLabel> <weight> ...
Option |
Type / Default |
Description |
|---|---|---|
recordsFile |
str, default=None |
Path (relative to the input file) to the flattened hanging-node records: one line per slave ‘<slaveLabel> <masterLabel> <weight> …’ with independent masters. Omit for dynamic AMR, where the adaptivity manager sets the records in memory via updateRecords(). |
testfiles/marmot/AMR_PatchTestU/test.inp*material, name=linearelastic, id=mat
30000, 0.15
*section, name=sec, material=mat, type=solid
all
*constraint, type=hangingnode, name=hn
recordsFile=hanging.dat
*job, name=amr, domain=3d
*solver, solver=NIST, name=S
*fieldOutput
>>perNode, elSet=all, field=displacement, result=U, name=U
*step, solver=S
maxInc=1.0, minInc=1e-4, maxNumInc=50, maxIter=25,
>>dirichlet, name=x0, nSet=x0, field=displacement, 1=0.0
>>dirichlet, name=x1, nSet=x1, field=displacement, 1=0.02
>>dirichlet, name=yf, nSet=yf, field=displacement, 2=0.0
>>dirichlet, name=zf, nSet=zf, field=displacement, 3=0.0
*Node
1, 0.0, 1.0, 1.0
2, 0.0, 0.0, 1.0
3, 0.0, 0.0, 0.0
4, 0.0, 1.0, 0.0
5, 1.0, 1.0, 1.0
6, 1.0, 0.0, 1.0
7, 1.0, 0.0, 0.0
8, 1.0, 1.0, 0.0
9, 0.0, 0.5, 1.0
10, 0.0, 0.0, 0.5
11, 0.0, 0.5, 0.0
12, 0.0, 1.0, 0.5
13, 1.0, 0.5, 1.0
14, 1.0, 0.0, 0.5
15, 1.0, 0.5, 0.0
16, 1.0, 1.0, 0.5
17, 0.5, 1.0, 1.0
18, 0.5, 0.0, 1.0
19, 0.5, 0.0, 0.0
20, 0.5, 1.0, 0.0
21, 2.0, 1.0, 1.0
22, 2.0, 0.0, 1.0
23, 2.0, 0.0, 0.0
24, 2.0, 1.0, 0.0
25, 2.0, 0.5, 1.0
26, 2.0, 0.0, 0.5
27, 2.0, 0.5, 0.0
28, 2.0, 1.0, 0.5
29, 1.5, 1.0, 1.0
30, 1.5, 0.0, 1.0
31, 1.5, 0.0, 0.0
32, 1.5, 1.0, 0.0
33, 1.0, 0.5, 0.5
34, 1.5, 0.5, 0.5
35, 1.5, 0.0, 0.5
36, 1.5, 0.5, 0.0
37, 1.0, 0.25, 0.5
38, 1.0, 0.0, 0.25
39, 1.0, 0.25, 0.0
40, 1.0, 0.5, 0.25
41, 1.5, 0.25, 0.5
42, 1.5, 0.0, 0.25
43, 1.5, 0.25, 0.0
44, 1.5, 0.5, 0.25
45, 1.25, 0.5, 0.5
46, 1.25, 0.0, 0.5
47, 1.25, 0.0, 0.0
48, 1.25, 0.5, 0.0
49, 1.5, 0.5, 1.0
50, 1.0, 0.25, 1.0
51, 1.0, 0.0, 0.75
52, 1.0, 0.5, 0.75
53, 1.5, 0.25, 1.0
54, 1.5, 0.0, 0.75
55, 1.5, 0.5, 0.75
56, 1.25, 0.5, 1.0
57, 1.25, 0.0, 1.0
58, 1.5, 1.0, 0.5
59, 1.0, 0.75, 0.5
60, 1.0, 0.75, 0.0
61, 1.0, 1.0, 0.25
62, 1.5, 0.75, 0.5
63, 1.5, 0.75, 0.0
64, 1.5, 1.0, 0.25
65, 1.25, 1.0, 0.5
66, 1.25, 1.0, 0.0
67, 1.0, 0.75, 1.0
68, 1.0, 1.0, 0.75
69, 1.5, 0.75, 1.0
70, 1.5, 1.0, 0.75
71, 1.25, 1.0, 1.0
72, 2.0, 0.5, 0.5
73, 2.0, 0.25, 0.5
74, 2.0, 0.0, 0.25
75, 2.0, 0.25, 0.0
76, 2.0, 0.5, 0.25
77, 1.75, 0.5, 0.5
78, 1.75, 0.0, 0.5
79, 1.75, 0.0, 0.0
80, 1.75, 0.5, 0.0
81, 2.0, 0.25, 1.0
82, 2.0, 0.0, 0.75
83, 2.0, 0.5, 0.75
84, 1.75, 0.5, 1.0
85, 1.75, 0.0, 1.0
86, 2.0, 0.75, 0.5
87, 2.0, 0.75, 0.0
88, 2.0, 1.0, 0.25
89, 1.75, 1.0, 0.5
90, 1.75, 1.0, 0.0
91, 2.0, 0.75, 1.0
92, 2.0, 1.0, 0.75
93, 1.75, 1.0, 1.0
*Element, type=C3D20
1, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20
2, 33, 14, 7, 15, 34, 35, 31, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48
3, 13, 6, 14, 33, 49, 30, 35, 34, 50, 51, 37, 52, 53, 54, 41, 55, 56, 57, 46, 45
4, 16, 33, 15, 8, 58, 34, 36, 32, 59, 40, 60, 61, 62, 44, 63, 64, 65, 45, 48, 66
5, 5, 13, 33, 16, 29, 49, 34, 58, 67, 52, 59, 68, 69, 55, 62, 70, 71, 56, 45, 65
6, 34, 35, 31, 36, 72, 26, 23, 27, 41, 42, 43, 44, 73, 74, 75, 76, 77, 78, 79, 80
7, 49, 30, 35, 34, 25, 22, 26, 72, 53, 54, 41, 55, 81, 82, 73, 83, 84, 85, 78, 77
8, 58, 34, 36, 32, 28, 72, 27, 24, 62, 44, 63, 64, 86, 76, 87, 88, 89, 77, 80, 90
9, 29, 49, 34, 58, 21, 25, 72, 28, 69, 55, 62, 70, 91, 83, 86, 92, 93, 84, 77, 89
*elSet, elSet=all, generate
1, 9, 1
*nSet, nSet=x0
1, 2, 3, 4, 9, 10, 11, 12
*nSet, nSet=x1
21, 22, 23, 24, 25, 26, 27, 28
72, 73, 74, 75, 76, 81, 82, 83
86, 87, 88, 91, 92
*nSet, nSet=yf
1, 2, 3, 4, 5, 6, 7, 8
10, 12, 14, 16, 17, 18, 19, 20
21, 22, 23, 24, 26, 28, 29, 30
31, 32, 35, 42, 46, 47, 54, 57
58, 64, 65, 66, 70, 71, 74, 78
79, 82, 85, 88, 89, 90, 92, 93
*nSet, nSet=zf
1, 2, 3, 4, 5, 6, 7, 8
9, 11, 13, 15, 17, 18, 19, 20
21, 22, 23, 24, 25, 27, 29, 30
31, 32, 36, 43, 47, 48, 49, 53
56, 57, 63, 66, 69, 71, 75, 79
80, 81, 84, 85, 87, 90, 91, 93
Implementing your own constraints
Subclass from the constraint base class in module edelweissfe.constraints.base.constraintbase
and implement its computeKernels. A constraint enters the equation system like an element does
(see The residual and the tangent): it adds its residual contribution \(\mathbf{P}\)
– the negative of the force it exerts on its nodes – and the tangent
\(\mathbf{K} = \partial \mathbf{P} / \partial \mathbf{U}\). Two examples:
A penalty spring of stiffness \(k\) that pulls a degree of freedom \(u\) towards \(\bar{u}\) exerts the force \(-k\,(u - \bar{u})\) on its node, so it adds \(P \mathrel{+}= k\,(u - \bar{u})\) and \(K \mathrel{+}= k\), exactly as a spring element would.
directionalspringpenaltyis this constraint.A Lagrange multiplier \(\boldsymbol{\lambda}\) enforcing \(\mathbf{g}(\mathbf{u}) = \mathbf{0}\) adds \(\boldsymbol{\lambda}^T \partial \mathbf{g} / \partial \mathbf{u}\) to the rows of the nodes and \(\mathbf{g}\) to the rows of the multipliers, with the tangent \(\begin{bmatrix} \boldsymbol{\lambda}^T \partial^2 \mathbf{g} / \partial \mathbf{u}^2 & (\partial \mathbf{g} / \partial \mathbf{u})^T \\ \partial \mathbf{g} / \partial \mathbf{u} & \mathbf{0} \end{bmatrix}\).
equalvaluelagrangianandrigidbodyare such constraints.
Because \(\mathbf{K}\) is the derivative of exactly what is added to \(\mathbf{P}\), a
constraint can be checked by finite differences of its own computeKernels: a wrong tangent, or a
wrong sign of \(\mathbf{P}\) alone, shows up as a mismatch. \(\mathbf{P}\) and \(\mathbf{K}\) negated
together pass this check, so the test of a constraint also checks the direction of the forces it exerts on its nodes:
a penalty pulls them back towards what it prescribes, a contact pushes the bodies apart, a Lagrange multiplier acts on
them as \(\boldsymbol{\lambda}^T \partial \mathbf{g} / \partial \mathbf{u}\) with the constraint’s own
\(\mathbf{g}\). tests/test_residual_tangent_convention.py does both for the constraints shipped with
EdelweissFE (the two integrated surface contact constraints in their own test modules); add a case there for a new
one. The helpers it uses are part of the package, so that the tests of a package with constraints of its own use them
as well:
Helpers for checking the rule every constraint follows: it adds a residual contribution P and the tangent
K = dP/dU (see the solver documentation, “The residual and the tangent”).
Because K must be the derivative of what is added to P, comparing it with central differences of P checks
both at once: a wrong tangent, or a P whose sign alone is wrong, shows up as a mismatch. A P and a K negated
together pass, so a test checks the direction of the forces the constraint exerts as well, which only the mechanics of
the constraint can say: that a penalty pulls its nodes back, that a contact pushes the bodies apart.
They are part of the package, so that the tests of a package that adds constraints of its own can check them the same way as the tests of EdelweissFE check its constraints.
- edelweissfe.utils.residualtangentcheck.centralDifferences(function, x, h=1e-07)[source]
The Jacobian of
functionatxby central differences.- Parameters:
function – Maps a vector to a vector.
x (ndarray) – Where to differentiate.
h (float) – The perturbation.
- Returns:
The matrix
d function / d x.- Return type:
np.ndarray
- edelweissfe.utils.residualtangentcheck.finiteDifferenceTangent(constraint, U, dU, timeStep=None, h=1e-07)[source]
The derivative of a constraint’s residual contribution by central differences.
UanddUare perturbed together: within an increment the unknown is the increment,U = U_n + dU, and a constraint with a history (friction) reads both.- Parameters:
constraint – The constraint.
U (ndarray) – The local solution.
dU (ndarray) – The local solution increment.
timeStep – The time step handed to the constraint.
h (float) – The perturbation.
- Returns:
The dense
nDof x nDofmatrixdP/dU.- Return type:
np.ndarray
- edelweissfe.utils.residualtangentcheck.localVector(constraint, nodalValue, scalarValues=())[source]
A constraint’s local degree-of-freedom vector, in the constraint’s own layout.
- Parameters:
constraint – The constraint.
nodalValue –
nodalValue(node, field)returns the values offieldatnode.scalarValues – The values of the constraint’s own scalar variables (e.g. Lagrange multipliers), which come last.
- Returns:
The local vector of size
constraint.nDof.- Return type:
np.ndarray
- Raises:
ValueError – If the values do not fill the constraint’s layout: their number is not
constraint.nDof.
- edelweissfe.utils.residualtangentcheck.residualContributionAndTangent(constraint, U, dU, timeStep=None)[source]
The residual contribution of a constraint and its tangent, the latter as one dense matrix.
The tangent is scattered from the flat VIJ values with the constraint’s own
(I, J)pattern – the one the global assembly uses – so that the check covers the pattern together with the values.- Parameters:
constraint – The constraint.
U (ndarray) – The local solution.
dU (ndarray) – The local solution increment.
timeStep – The time step handed to the constraint.
- Returns:
Pand the densenDof x nDoftangent.- Return type:
tuple[np.ndarray, np.ndarray]
A constraint takes part in the equation system while it is active. A constraint that is not active gets no
degrees of freedom, the solvers do not evaluate it, and its scalar variables – the Lagrange multipliers of
equalvaluelagrangian, say – leave the equation system with it, so that they leave no unknowns without an equation
behind. Two states decide it, and active is true while both are:
switchedOn: the constraint is switched on from outside. Settingactiveswitches it, e.g. amodelupdatestep action that switches a constraint off for a step, as intestfiles/marmot/PenaltyBasedIndirectControl. A switch-off from outside lasts until the constraint is switched on again from outside.withinReach: the constraint finds itself within reach of what it constrains, in itsupdateConnectivity. A contact constraint whose partner is out of reach sets it to false, and reports the change ofactiveas a change of its connectivity, so that the solver rebuilds its equation system. A constraint that does not check stays within reach.
So a constraint switched off from outside stays off, whatever it finds in reach; switched on again, it takes part
while it is within reach. Switched off from outside, a constraint rests: the model does not locate it anew in
updateConnectivity, which it does for every constraint that is switched on, those out of reach included, so that
they can find themselves within reach again. A constraint that is not active accepts no state either
(acceptLastState): a contact switched off for a step keeps its multipliers and its frictional history as they
were. A solver that cannot handle a constraint, e.g. NEST a contact, refuses it only while it is switched on.
- class edelweissfe.constraints.base.constraintbase.ConstraintBase(name, model, *args, **kwargs)[source]
The constraint base class.
Constraints can act on nodal variables, and scalar variables. If scalar variables are required, the can be created on demand by defining
getNumberOfAdditionalNeededScalarVariables()andassignAdditionalScalarVariables(), which are called at the beginning of an analysis.If scalar variables are used, EdelweissFE expects the layout of the residual contribution
P(and of the tangent) to be of the form[ node 1 - dofs field 1, node 1 - dofs field 2, node 1 - ... , node 1 - dofs field n, node 2 - dofs field 1, ... , node N - dofs field n, scalar variable 1, scalar variable 2, ... , scalar variable J ].
- Parameters:
name (str) – The name of the constraint.
model (FEModel) – A dictionary containing the model tree.
kwargs – Key value pairs.
- Attributes:
activeWhether this constraint takes part in the equation system: while it is
switchedOnandwithinReach.- checkpointedState
couplesAllItsDofsWhether the entity’s block couples every one of its degrees of freedom with every other one.
fieldsOnNodesThe fields on the nodes this constraint is acting on.
nDofThe total number of degrees of freedom this constraint is associated with.
nodesThe nodes this constraint is acting on.
- schema
switchedOnWhether this constraint is switched on from outside, e.g.
Methods
Called by
advanceToTime()when an increment is accepted, so a stateful constraint (e.g. frictional contact) can promote the state of the last (converged) Newton iterate to its history.assignAdditionalScalarVariables(scalarVariables)Assign a list of scalar variables associated with this constraint.
computeKernels(U_np, dU, P, K, timeStep)Add this constraint's contribution
Pto the residual, and its tangentK = dP/dU.computeKernelsExplicit(U_np, dU, P, timeStep)Add this constraint's contribution to the residual, without a tangent.
fromConstraintDefinition(name, definition, model)Create this constraint from a parsed
.inpconstraint definition.This method is called to determine the additional number of scalar variables this Constraint needs.
The state this constraint carries from one increment to the next; by default the attributes declared in
checkpointedState.getVIJContributionSize()Return the number of entries this entity contributes to the VIJ (COO) system matrix.
initializeVIJContribution(idcs, I_, J_, offset)Fill the global
IandJindex arrays for this constraint's VIJ contribution.setRestartData(data)Restore the state
getRestartData()returned.shapeVIJContribution(flat_view)Shape the flat VIJ values slice for this entity.
updateConnectivity(model)Called once at the start of every increment, before the global equation system is (re)created, so a constraint can refresh a dynamic candidate set (e.g. proximity-based contact pairs) and declare whether its
nodes/fieldsOnNodes/nDoffootprint has changed since the last call.- classmethod fromConstraintDefinition(name, definition, model, journal=None)[source]
Create this constraint from a parsed
.inpconstraint definition.The one place a module’s input-file shape (set/surface names, string-typed booleans) is turned into the typed arguments its real constructor takes. Override it together with a typed
__init__; leave it alone and the dict-consuming constructor is used unchanged.- Parameters:
- Returns:
The constructed constraint.
- Return type:
- acceptLastState()[source]
Called by
advanceToTime()when an increment is accepted, so a stateful constraint (e.g. frictional contact) can promote the state of the last (converged) Newton iterate to its history. Only anactiveconstraint is called: one that took no part in the increment keeps its history.The default implementation does nothing, which is correct for every stateless constraint (i.e. every constraint that does not override this method).
- assignAdditionalScalarVariables(scalarVariables)[source]
Assign a list of scalar variables associated with this constraint.
- Parameters:
list – The list of ScalarVariables associated with this constraint.
scalarVariables (list[ScalarVariable])
- abstractmethod computeKernels(U_np, dU, P, K, timeStep)[source]
Add this constraint’s contribution
Pto the residual, and its tangentK = dP/dU.A constraint follows the rule every element and every load follows (see The residual and the tangent):
Pis the force the constraint needs at its degrees of freedom, i.e. the negative of the force it exerts on the nodes, exactly as an element of the same mechanics would add its internal force. A penalty spring of stiffness \(k\) that pulls a degree of freedom \(u\) towards \(\bar{u}\), for example, adds\[P \mathrel{+}= k\,(u - \bar{u}), \qquad K \mathrel{+}= k .\]A Lagrange multiplier constraint \(g(\boldsymbol{u}) = 0\) adds \(\boldsymbol{\lambda}^T \partial g / \partial \boldsymbol{u}\) to the rows of its nodes and \(g\) to the rows of its multipliers.
Both
PandKare augmented, never overwritten.Kmust be the derivative of what is added toP: a test can then check it by finite differences, which also catches aPwhose sign alone is wrong. APand aKnegated together pass that check, so a test checks the direction of the forces the constraint exerts as well.- Parameters:
U_np (ndarray) – The current solution, restricted to this constraint’s degrees of freedom; its own scalar variables (e.g. Lagrange multipliers) come last.
dU (ndarray) – The solution increment since the last converged state, likewise restricted.
P (ndarray) – The residual contribution to augment.
K (ndarray) – The tangent
dP/dUto augment, in the shapeshapeVIJContribution()gives it (a densenDof x nDofview unless the constraint declares a sparser pattern). A constraint that inheritsForcesOnlyExplicitEvaluationis also called withK=Noneand must then skip the tangent.timeStep (TimeStep) – The current time step.
- computeKernelsExplicit(U_np, dU, P, timeStep)[source]
Add this constraint’s contribution to the residual, without a tangent.
The counterpart of
computeKernels()for an explicit solver, which solves no linear system and so has no use for a tangent. This default callscomputeKernels()with a scratch tangent that is thrown away, so that every constraint also works in an explicit analysis. A constraint whose tangent is costly overrides it, e.g. throughForcesOnlyExplicitEvaluation.- Parameters:
U_np (ndarray) – The current solution, restricted to this constraint’s degrees of freedom.
dU (ndarray) – The current solution increment, likewise restricted.
P (ndarray) – The residual contribution to augment, with the sign described in
computeKernels().timeStep (TimeStep) – The current time step.
- getNumberOfAdditionalNeededScalarVariables()[source]
This method is called to determine the additional number of scalar variables this Constraint needs.
- Returns:
Number of ScalarVariables required.
- Return type:
int
- getRestartData()[source]
The state this constraint carries from one increment to the next; by default the attributes declared in
checkpointedState. Overridden where the state is not a plain attribute.- Returns:
A flat mapping of array name to array.
- Return type:
dict[str, numpy.ndarray]
- setRestartData(data)[source]
Restore the state
getRestartData()returned.- Parameters:
data (dict[str, ndarray]) – The mapping of arrays.
- updateConnectivity(model)[source]
Called once at the start of every increment, before the global equation system is (re)created, so a constraint can refresh a dynamic candidate set (e.g. proximity-based contact pairs) and declare whether its
nodes/fieldsOnNodes/nDoffootprint has changed since the last call.The default implementation does nothing and reports no change, which is correct for every constraint whose DOF footprint is fixed at construction (i.e. every constraint that does not override this method).
- Parameters:
model (FEModel) – The current model.
- Returns:
True if this constraint’s contribution to the global system has changed and the equation system must be rebuilt before the next Newton solve.
- Return type:
bool
- property active: bool
while it is
switchedOnandwithinReach. A constraint that does not take part gets no degrees of freedom and adds nothing to the residual or the tangent; its scalar variables (e.g. Lagrange multipliers) leave the equation system with it, so that they leave no empty rows behind.Setting it switches the constraint on or off from outside, e.g. by a
modelupdatestep action at the start of a step. A switch-off from outside lasts until the constraint is switched on again from outside, whatever the constraint finds in reach; switched on again, it takes part while it is within reach.- Returns:
Whether the constraint takes part in the equation system.
- Return type:
bool
- Type:
Whether this constraint takes part in the equation system
- checkpointedState: dict | None = None
The state this constraint carries from one increment to the next, by attribute name and type; see
checkpointedstate. Every constraint declares it (empty if it carries none) or overridesgetRestartData()andsetRestartData().
- abstract property fieldsOnNodes: list
The fields on the nodes this constraint is acting on.
- abstract property nDof: int
The total number of degrees of freedom this constraint is associated with.
- abstract property nodes: list
The nodes this constraint is acting on.
- property switchedOn: bool
Whether this constraint is switched on from outside, e.g. by a
modelupdatestep action, which switches it throughactive. A constraint switched off from outside rests: the model does not locate it anew (updateConnectivity()), it accepts no state (acceptLastState()), so that its history stays as it was, and the solvers neither evaluate it nor refuse it. A solver that cannot handle a constraint refuses it while it is switched on, even out of reach: it may come within reach during the step.- Returns:
Whether the constraint is switched on.
- Return type:
bool
- withinReach: bool = True
Whether this constraint is within reach of what it constrains, as it finds itself in
updateConnectivity(): e.g. a contact constraint whose partner is far away is not, and takes no part in the equation system then (seeactive). A constraint sets it there, and reports a change ofactiveas a change of its connectivity, so that the solver rebuilds its equation system. A constraint that does not check stays within reach.
Multi-point (DOF-elimination) constraints subclass from
edelweissfe.constraints.base.multipointconstraintbase instead – they contribute nothing to
the load vector or system matrix; they only declare linear dependency records, which the solvers
condense via MultiPointConstraintTransformation.
- class edelweissfe.constraints.base.multipointconstraintbase.MultiPointConstraintBase(name, model, *args, **kwargs)[source]
The multi-point constraint base class.
- Parameters:
name (str) – The name of the constraint.
model (FEModel) – The model tree.
kwargs – Key value pairs of options.
- Attributes:
fieldsOnNodesThe fields this constraint needs on each of
nodes, as forfieldsOnNodes.nodesThe nodes this constraint needs fields on, as for
nodes.- schema
Methods
Called by
advanceToTime()when an increment is accepted, so a stateful multi-point constraint can promote the state of the last (converged) iterate to its history.Return the nodes this constraint claims as multi-point-constraint slaves.
fromConstraintDefinition(name, definition, model)Create this constraint from a parsed
.inpconstraint definition.getMultiPointConstraints(dofManager)Return the linear dependency records of this constraint in global degree-of-freedom indices of the given DofManager.
- classmethod fromConstraintDefinition(name, definition, model, journal=None)[source]
Create this constraint from a parsed
.inpconstraint definition.See
edelweissfe.constraints.base.constraintbase.ConstraintBase.fromConstraintDefinition()for the rationale; this is the same seam for the multi-point-constraint hierarchy.- Parameters:
- Returns:
The constructed constraint.
- Return type:
- acceptLastState()[source]
Called by
advanceToTime()when an increment is accepted, so a stateful multi-point constraint can promote the state of the last (converged) iterate to its history.The default implementation does nothing, which is correct for every stateless constraint (i.e. every constraint that does not override this method).
- claimedSlaveNodes()[source]
Return the nodes this constraint claims as multi-point-constraint slaves.
Peer constraints use this to avoid claiming the same slave twice, which the condensation operator would reject. The default derives the set from
_records, whose entries are expected to carry the slave node as their first item.- Returns:
The slave nodes of this constraint.
- Return type:
set
- abstractmethod getMultiPointConstraints(dofManager)[source]
Return the linear dependency records of this constraint in global degree-of-freedom indices of the given DofManager.
- Parameters:
dofManager – The current
DofManager, used to resolve (node, field, component) to global degree-of-freedom indices.- Returns:
One record per slave degree of freedom:
(slaveDofIndex, [(masterDofIndex, coefficient), ...]).- Return type:
list[tuple[int, list[tuple[int, float]]]]
- property fieldsOnNodes: list
The fields this constraint needs on each of
nodes, as forfieldsOnNodes.
- mayYieldSlaveToEarlierClaim: bool = True
Whether a slave DOF of this constraint may be left to an earlier constraint that claims it too (its record is then dropped). If False (e.g. hanging nodes), a dropped record is an error.