Unable to receive traction from two different coupled surfaces

In the following scenario, we have a left muscle that is activated and contracts. This should result in the tendon moving to the left and that the right muscle (a passive muscle) elongating. The transfer of information is such, that the tendon reads the traction sent by both the left and right muscle.

muscle_tendon_muscle

Contrary to what is expected, the results of the simulation show that the tendon and the right muscle do not move. The tendon receives non-zero traction values from the left muscle, but the internal traction field of the muscle remains zero. Note that in the plot we only show the z coordinate of materialTraction.

option3_drama1

The terminal output of the tendon participant shows why the traction field remains zero (I added some comments for a better understanding):

**dynamic hyperelasticity updates boundary condition**

**input from the left muscle**

PreCICE surface coupling, timestep 61, t=0.61
read and set Neumann BC: 
{el. 0, "2-", dofVectors: [(0: (1.40108e-06,1.97848e-06,-3.84902e-06)),(1: (2.46872e-07,1.68348e-06,-3.90705e-06)),(2: (-1.73185e-06,2.1483e-06,-4.57642e-06)),(3: (1.77552e-06,1.022e-07,-2.90253e-06)),(4: (1.23194e-07,-3.40221e-07,-1.82769e-06)),(5: (-1.76219e-06,-4.38139e-07,-3.05034e-06)),(6: (2.61164e-06,-2.21063e-06,-6.22083e-06)),(7: (-2.93468e-07,-1.97306e-06,-2.11356e-06)),(8: (-2.09342e-06,-1.53145e-06,-5.33528e-06)),], surfaceDofs: 0,1,2,3,4,5,6,7,8,} 

update neumann BC in hyperelasticity solver 
component 0, neumannBoundaryConditions_ rhs values: [3.75879e-07,1.85336e-06,5.82036e-07,2.08081e-07,3.12227e-07,-1.39291e-07,-3.69225e-07,-1.7606e-06,-5.56101e-07,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
component 1, neumannBoundaryConditions_ rhs values: [4.81205e-07,-1.5853e-08,-5.38283e-07,1.67579e-06,-1.05575e-06,-2.03633e-06,4.27177e-07,-3.63224e-07,-4.64973e-07,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
component 2, neumannBoundaryConditions_ rhs values: [-7.51902e-07,-2.8026e-06,-1.00897e-06,-3.42823e-06,-9.30967e-06,-2.30651e-06,-9.56441e-07,-2.90504e-06,-8.10272e-07,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]

**input from the right muscle**

read and set Neumann BC: 
{el. 1, "2+", dofVectors: [(36: (-1.05181e-30,-1.05181e-30,8.88178e-16)),(37: (0,-1.05181e-30,8.88178e-16)),(38: (-1.05181e-30,-1.05181e-30,8.88178e-16)),(39: (-1.05181e-30,0,8.88178e-16)),(40: (0,0,8.88178e-16)),(41: (-1.05181e-30,0,8.88178e-16)),(42: (-1.05181e-30,-1.05181e-30,8.88178e-16)),(43: (0,-1.05181e-30,8.88178e-16)),(44: (-1.05181e-30,-1.05181e-30,8.88178e-16)),], surfaceDofs: 18,19,20,21,22,23,24,25,26,} 

update neumann BC in hyperelasticity solver 
component 0, neumannBoundaryConditions_ rhs values: [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1.57772e-31,-2.10363e-31,-1.57772e-31,-6.31089e-31,-8.41452e-31,-6.31089e-31,-1.57772e-31,-2.10363e-31,-1.57772e-31]
component 1, neumannBoundaryConditions_ rhs values: [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1.57772e-31,-6.31089e-31,-1.57772e-31,-2.10363e-31,-8.41452e-31,-2.10363e-31,-1.57772e-31,-6.31089e-31,-1.57772e-31]
component 2, neumannBoundaryConditions_ rhs values: [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2.22045e-16,8.88178e-16,2.22045e-16,8.88178e-16,3.55271e-15,8.88178e-16,2.22045e-16,8.88178e-16,2.22045e-16]
DynamicHyperelasticitySolver, timestep 0/1, t=0.61

**computation of the dynamic hyperelasticity solver starts**

...

The boundary condition from the right solver overwrites the boundary condition from the left solver (see final rhs values).

As this point, we ask us the following:

Can we avoid the situation where a tendon receives traction from two different coupled surfaces?

Yes, as long it is possible to send the traction from the tendon to the muscle. For example:

muscle_tendon_muscle(1)

At the moment this is not yet trully possible (see #18) .

How can we change the code in order to support receiving traction from two different boundary conditions?

It may help to think how of how we do that if we have a single tendon where we manually provide boundary conditions at different faces. Then we also have to "append" the boundary conditions:

# neumann bc
k = mz-1

variables.elasticity_neumann_bc = [{"element": j*mx + i, "constantVector": [0,0, 5], "face": "2-", "isInReferenceConfiguration": True} for j in range(my) for i in range(mx)]

 for j in range(my):
   for i in range(mx):
     variables.elasticity_neumann_bc.append({"element": k*mx*my + j*mx + i, "constantVector": [0,0, 1], "face": "2+", "isInReferenceConfiguration": True})

Why is it that the tendon always reads the boundary condition from the left first? I have tried different ordering in different coupling schemes: multi-coupling or two explicit schemes. Also I have tried using a dynamic hyperelasticity solver instead of the MuscleContractionSolver, but I didn't manage to change the order in which the tendon reads the tractions.

maierbn / opendihu

Unable to receive traction from two different coupled surfaces #17