In this commit, the traction element was created to set the traction force to zero at the nodes that are in contact within a specified contact interaction. This feature is required in the coupling procedure with Metalub (see commit 3186 of Metafor and commit 1456 of Metalub for more details).

A ContactTraction2DElement or ContactTraction3DElement is created like a Traction2DElement or Traction3DElement. The only additional property that has to be set is NBR_CONTACT_INTERACTION, which is the number of the contact interaction. For instance,

prp = ElementProperties(ContactTraction2DElement)
prp.put(PRESSURE, -p['pMax'])
prp.put(NBR_CONTACT_INTERACTION, 3)

Hence, if nodes in this ContactTraction2DElement are slave nodes in contact in the contact interaction 3, the traction force, which is applied to these nodes, is set equal to zero.

Instead of implementing the new feature directly in the TractionElement of the mtElement package, a new template class ContactTractionElement was created in the mtContact (elements) package. It had to be moved to the mtContact package (instead of the mtElement package) in order to have access to the contact features via include “ContactInteraction.h”. Moreover, the ContactTractionElement derives of TractionElement to minimize duplicate code. In a similar way to the TractionElement, 2D and 3D elements were created: ContactTraction2DElement and ContactTraction3DElement. This required to export all methods of Traction2DElement and Traction3DElement via MTELEMENTS_API (+ ELEMENT_CPP_MACRO2 extension), which only had to be added in Traction3DElement since it was already done for Traction2DElement, certainly to derive it in mtXFEM.

Some minor modifications:

• The new property NBR_CONTACT_INTERACTION was added in mtElements.h/cpp as an integer property.
• It was necessary to add ContactTraction2DElement and ContactTraction3DElement to the functions is2D and is3D in LoadingInteraction.cpp in order to be able to add them to the LoadInteraction in the Python data set.
• Untabify in ContactElement.h/cpp.
• getSlaveNode(), previously private method in ContactElement.h/cpp, becomes public to use it in ContactTractionElement.inl.
• The method isInContactPrev() was added in ContactElement.h/cpp to get the contact status at the end of the previous time iteration. Using this contact status has the advantage that the contact traction activation/deactivation does not happen during the same time step, i.e. during the mechanical iteration. It might be possible to use the current contact status but this could introduce some undesired oscillations.

Concerning the implementation, the ContactTractionElement is very similar to the TractionElement except from parts of the following virtual functions:

• checkProps: This function checks whether the new property NBR_CONTACT_INTERACTION was set during the initialization of the elements. If this property was set, it checks whether the specified number corresponds to an existing contact interaction. Once this interaction is found, pointers to the slave nodes corresponding to the nodes of the ContactionTractionElement are saved in the new ces private vector of Contact ElementS (pointer) in the ContactTractionElement.

• fillExternalForces: This function uses the fillExternalForces function of TractionElement first and then it sets all contact forces corresponding to nodes in contact equal to zero. As mentionned earlier, the contact state is determined based on the results of the previous time integration to reduce possible oscillations.

• noAnalyticalStiffnessAvailable: An analytical stiffness matrix was defined for the ContactTractionElement. Hence, this function can return false.

• fillAnalyticalStiffness: Similar to the fillExternalForces function, the fillAnalyticalStiffness of TractionElement is used first and then the rows corresponding to the traction forces that were set equal to zero are also set equal to zero. The structure of this matrix is illustrated in a comment of the fillAnalyticalStiffness function to understand the details.

Two test cases were created in the folder mtContact/tests to test the ContactTractionElement in 2D and 3D. Hence, they have the same names than the elements that are tested: contactTraction2DElement.py and contactTraction3DElement. The 3D test is the natural extension of the 2D case. Hence, only some details about the 2D case are explained here.

As shown in the previous figure, two squares are crushed by a compression pressure and a contact tool that moves downwards by $2 \; dy$. In the case on the left, the traditional Traction2DElement is used, while the new ContactTraction2DElement is used on the right. The ContactTraction2DElement works correctly, if the total force on node $P_{14}$ along $\mathbf{y}$ is equal to the contact force along $\mathbf{y}$ at the end of the computation, i.e. that the traction force has been set to zero for the node in contact. To check this equality, the difference between both forces (total and contact) is computation in the extractor diffTotalAndContactForceInCase2. The results of the cases with the Traction2DElement and the ContactTraction2DElement are also compared by computing the differences of the total force along $\mathbf{y}$ at the nodes $P_{4}$ and $P_{14}$ in the extractor diffTotalForcesBetween2Cases. This difference should be different from 0, unlike the previous one.