Implement multipoint constraint for EP simulations

[javijv4]

Problem

Description

For EP simulations that couple 3D and 1D simulations, the current approach is to use a projection that uses <End_nodes_face_file_path>. This works by assigning the closest 3D degree of freedom to the end node of the 1D mesh, effectively changing the problem's geometry. For example, the first image below shows the original 1D mesh (images taken from #474) and the second shows the mesh that svMultiPhysics is using to solve the problem (white is the 1D Purkinje mesh).

Reproduction Steps

See the MWE mwe.zip or the test in tests/cases/cep/cylinder_purkinje_1d3d

Expected Behavior

Ideally, we would solve without the need to change the 1D mesh, making the simulations more robust to changes in mesh resolution. One approach to do so is using linear Multi-Point Constraints (MPC), see here.

Solution

MPC works by constraining the DOFs of the 3D element that contains the 1D node to match the value at the position where the 1D node is. So, for a 3D element that contains a 1D node, we want,

$$\sum_{a=1}^{n_{nodes}^{3D}} N_a(x_{1D})^{3D} \phi_a^{3D} = \phi^{1D}(x_{1D})$$

where $N_a$ are the 3D surface basis functions, $\phi_a^{3D}$ are the nodal values of the 3D element, and $\phi^{1D}$ are the nodal values of the node to be constrained.

MPC matrix implementation

For an uncoupled 3D-1D simulation, at each time step, and at each Newton iteration, svMultiphysics solves,

We can augment that system to introduce the constraint through a Lagrange multiplier,

To avoid modifying the system matrix in the main solve loop, the MPC constraints can be imposed after the system matrix is assembled by solving for $\lambda$,

and modifying the RHS,

Implementation

To implement this, we'll need to

1. Create a function to read the MPC nodes and find the 3D element that contains each node.
2. Create a function that, after the regular assembly, modifies the RHS based on the equations above.
3. Deal with parallelization issues.

Additional context

No response

Code of Conduct

• I agree to follow this project's Code of Conduct and Contributing Guidelines

[ktbolt]

@javijv4 If I understand correctly you are planning to constrain the 3D elements DOF so that the nodes of the 1D Purkinje mesh don't move in time (i.e., cardiac cycle) ?

[javijv4]

@ktbolt No, but it would be cool to do something like that at some point.

I want to implement it for electrophysiology and heat physics, which are scalar equations that do not involve mesh movement. The idea is to constrain the 3D mesh at the position of a 1D constrained node to match the 1D value. In EP, the excitation starts at the 1D mesh, and the potential propagates through it until it reaches the end-nodes that need to be coupled to the 3D mesh.

Does that make sense?

[ktbolt]

@javijv4 Thank you for the clarification !