Labeling overlapping nodes automatically

[mallickrishg]

@brendanjmeade This is a problem related to automating the labeling of nodes of fault segments that intersect. Consider a geometry of a wavy fault in a rectangular box (attached image), I was wondering if you had suggestions for how to come up with a general scheme to label which fault nodes overlap between which fault elements? So far, this has been easy because I only dealt with sequential elements and the first and last fault elements had free/unconstrained nodes.

With closed loops, it's a bit different but possibly simpler since there are no unconstrained nodes - the first and last element end up sharing a node.

For intersecting faults however, I think this might be a problem. In the figure below, elements 32,33 & 38 have a common node and similarly 14,15 & 46. I will think about this today.

[brendanjmeade]

This is a problem related to automating the labeling of nodes of fault segments that intersect. Consider a geometry of a wavy fault in a rectangular box (attached image), I was wondering if you had suggestions for how to come up with a general scheme to label which fault nodes overlap between which fault elements? So far, this has been easy because I only dealt with sequential elements and the first and last fault elements had free/unconstrained nodes.

My initial response would be to focus on identifying (and matching) "hanging endpoints".

loop over each node that does not already have two matched endpoints by construction
    loop over each endpoint
        find inidices of elements with matching (same xy coordinates) endpoints (build big list of these)

It's $O(n^2)$ but should be fine for all but huge problems. Is this along the lines that you were thinking of?

With closed loops, it's a bit different but possibly simpler since there are no unconstrained nodes - the first and last element end up sharing a node.

For intersecting faults however, I think this might be a problem. In the figure below, elements 32,33 & 38 have a common node and similarly 14,15 & 46. I will think about this today.

This case is really interesting. I'd guess we'll have to learn how to write continuity of slip and slip gradient for three faults rather than two. Does that sound about right to you? Do you have a target notebook with this model geometry?

[mallickrishg]

Here is the notebook i am working on for this problem https://github.com/brendanjmeade/bemcs/blob/main/examples/examples_3_intersectingfaults.ipynb

This case is really interesting. I'd guess we'll have to learn how to write continuity of slip and slip gradient for three faults rather than two. Does that sound about right to you? Do you have a target notebook with this model geometry?

I was thinking that by writing slip continuity and smoothness equations 2 nodes at a time, we will be ok. So in the case of 3 fault segments meeting at a node, we write 2 sets of equations:

1. For segment 1 and 2 that meet at the common node
2. For segment 2 and 3 that meet at the same node That should automatically force segments 1 and 3 to close. The difficulty right now is in labelling those 'overlapped' nodes - sort of how triangulation of a set of vertices to produce triangles. I will try what you suggested above and keep you posted.

[mallickrishg]

@brendanjmeade I just realized something quite interesting. It's related to solving the above mentioned problem.

• If we have an open loop meshed with N elements, then there are N-1 internal overlapping nodes (need boundary conditions at the end node)
• If we have a simply closed loop with N elements, then there are N internal nodes (no end node boundary conditions needed)
• If we have 2 open loops that intersect at a point, with N and M elements respectively, then we have N+M-1 equations for internal nodes
• If we have 2 intersecting closed loops like in the figure at the top of this issue, then we end up with N+1 equations of overlapping nodes. There is a redundancy somewhere.

I need to think more carefully about where to apply continuity of slip and smoothness of slip conditions. Wow! this is tricky!

[brendanjmeade]

I've got to things about this more too. The way currently thinking about the problem that all we have to do is determine which element nodes touch each other and then we apply your boundary continuity relationships to the appropriate rows. In this perspective, I don't really have the concept of a loop. It's still there geometrically, but not algebraically. All there is are intersections. Does this make some sense?

[mallickrishg]

You're right, the idea is really centered on intersecting line segments

• a hanging/open node has only 1 element and requires extra equations to constrain slip
• typical interior nodes have 2 elements and do not require any extra equations
• triple junction nodes have 3 elements and end up overdetermined, I need to figure out if there is a redundancy or the extra equations are necessary. Will create a new example notebook for a 3 element triple junction and investigate this. This stuff takes me back to math olympiad camp.... it used to happen over summer too. So fun!

[mallickrishg]

Very interesting results. Turns out our history of working on block models really came to the rescue. For a triple junction we need to conserve the slip vector and the slip-gradient vector of the 3 segments evaluated at the common node in a clockwise (or anti-clockwise) loop just like we would for block modeling. I have implemented this in a new notebook - https://github.com/brendanjmeade/bemcs/blob/main/examples/examples_3_simpleintersectingfaults.ipynb One of the unexpected results of these calculations is now I have a system with fewer equations than unknowns because the overlapping node is a member of 3 elements and so results in 2 fewer equations than usual. I have not yet figured out what this means or if this is something that we need to fix but I am getting away with the matrix inverse by using pseudo inverse.

The resulting stress distributions are non singular - however they need not be continuous (see attached figure). This is reasonable since if you recall that $\sigma_{xx}$ for a horizontal fault is discontinuous but not singular on either side of the fault when there exists any slip gradient.

[brendanjmeade]

Wow! Now that's a figure of an interesting problem. I've read through the code a little bit to look at the matrix inverse part. Here's how I've reasoned about it.

It looks like the matrix is square but rank degenerate:

matrix_system.shape=(18, 18)
np.linalg.matrix_rank(matrix_system)=16
np.linalg.cond(matrix_system)=1.7796848236864911e+18

It should be square (as all BEM problems are!). The question is then, "How is the matrix rank degenerate?". Some thoughts:

• Looking at matrix_system it looks like some of the rows are all zeros. Specifically, rows 7 and 12. Should these be all zeros? I'm sort of guessing no because it would mean that they do not contribute to the solution in any way. This also means that the linear system is essentially missing two equations which might account for the difference between the matrix size (18) and the matrix rank (16)
• Also, rows 6 and 13 each sum almost exactly to zero. Not sure if this is by design or not?

Original matrix_system

Row normalized matrix_system

[mallickrishg]

• Looking at matrix_system it looks like some of the rows are all zeros. Specifically, rows 7 and 12. Should these be all zeros? I'm sort of guessing no because it would mean that they do not contribute to the solution in any way. This also means that the linear system is essentially missing two equations which might account for the difference between the matrix size (18) and the matrix rank (16)

Exactly. Those two rows of zeros is by construction because I expected there to be 2 more equations than I was able to construct.

• Also, rows 6 and 13 each sum almost exactly to zero. Not sure if this is by design or not?

These are the rows where I do the triple junction kinematics thing. s_1 - s_2 - s_3 = 0 and similarly for the derivative of slip.

Putting it all together, it is strange that we end up with 2 fewer equations per triple junction than other cases, but the system is still closed and we end up with a unique solution. I haven't tested enough cases to be convinced that it is a unique solution but it is a solution. I've added some description about the number of equations in the design matrix and how it changes with number of triple junctions here - https://github.com/brendanjmeade/bemcs/blob/main/examples/examples_3_simpleintersectingfaults.ipynb

[brendanjmeade]

I think I get it now. There are only 16 boundary conditions and 18 quadratic coefficients to solve for. That means that the problem is underdetermined. You inserted the blank rows to keep the matrix square.

Exactly. Those two rows of zeros is by construction because I expected there to be 2 more equations than I was able to construct.

I so wish we had a whiteboard and a few hours right now! Three total guesses:

1. Are there two terms in the state vector (quadratic_coefs) that should be identical by construction? If so, we can reduce the number of columns by two.
2. Do we need to triple junction kinematics equations? Is slip continuity alone sufficient to ensure no singularities?
3. Alternatively, could we write the kinematics equations as more than one equation?

\begin{bmatrix}
   0 \\
   0
\end{bmatrix}
= 
\begin{bmatrix}
   1 & -\frac{1}{2} & 0\\
   0 & -\frac{1}{2} & -1
\end{bmatrix}
=
\begin{bmatrix}
   s_1 \\
   s_2 \\
   s_3
\end{bmatrix}

instead of,

   0 
= 
\begin{bmatrix}
   1 & -1 & -1
\end{bmatrix}
=
\begin{bmatrix}
   s_1 \\
   s_2 \\
   s_3
\end{bmatrix}

I don't think that's right.

[mallickrishg]

I don't think breaking it up into 2 equations will work. We really need to sit together and whiteboard this. Are you available tomorrow?

[brendanjmeade]

I don't think breaking it up into 2 equations will work.

I agree that this is a bad idea.

We really need to sit together and whiteboard this. Are you available tomorrow?

Yes! After 12 EST. Is there a time that works for you?

[mallickrishg]

That works for me. I'm free all day. 12.05pm EST is good for me. Same zoom link?

[brendanjmeade]

Reading through the matrix assembly and trying a quick experiment. What about replacing,

matrix_system[4, :] = matrix_slip[4, :] - matrix_slip[6, :] - matrix_slip[12, :]  # x component
matrix_system[5, :] = matrix_slip[5, :] - matrix_slip[7, :] - matrix_slip[13, :]  # y component

with something like,

# elements 1 and 2
matrix_system[4, :] = matrix_slip[4, :] - matrix_slip[6, :]  # x component
matrix_system[5, :] = matrix_slip[5, :] - matrix_slip[7, :]  # y component
# elements 1 and 3
matrix_system[6, :] = matrix_slip[6, :] - matrix_slip[12, :]  # x component
matrix_system[7, :] = matrix_slip[7, :] - matrix_slip[13, :]  # y component
# elements 2 and 3
matrix_system[12, :] = matrix_slip[12, :] - matrix_slip[4, :]  # x component
matrix_system[13, :] = matrix_slip[13, :] - matrix_slip[5, :]  # y component

PRETTY SURE I GOT ALL OF THE INDICES WRONG!!!

And similar for the slip gradient? It'll definitely make things full rank.

[mallickrishg]

That is a good idea and it was what I tried to get to work all day on Friday but there are 2 problems that come up -

1. With only continuity of slip, stress singularities were still there. Although now I can't recall implementing all 3 sets of equations. Something for me to check tomorrow morning.
2. If we implement only 4 of those 6 equations (only 4 are needed because if a = b, b = c then a = c) for slip and similarly for slip gradient, we end up with 2 more equations and so now the system is overdetermined. In that case the pseudo-inverse solution ends up with substantial residuals. This suggests there is no redundancy in these equations and they conflict with each other
3. Based on the stress calculations for

   matrix_system[4, :] = matrix_slip[4, :] - matrix_slip[6, :] - matrix_slip[12, :]  # x component
   matrix_system[5, :] = matrix_slip[5, :] - matrix_slip[7, :] - matrix_slip[13, :]  # y component

the singularities disappear even though slip at a given triple junction node is not exactly the same value for all 3 segments. But I am finding some issues regarding how slip is partitioned among the 3 fault elements for a limiting case where I set one of the faults to have 0 slip, the trivial solution would be to just treat it as if that fault does not exist and it is a simple 2-fault system. But because we are missing a couple equations, the solution does not return this answer.

[brendanjmeade]

This turns out to be quite the geometry/kinematics/bookkeeping/linear algebra problem. That feels like the heart of most problems in this field. I'd like to understand your point 2 above better when we chat. Also, one random and unprincipled idea: What happens if we modify the model a little bit? Instead of just the three faults, what if we put the three faults inside a triangle that bounds the three interior elements? This would create 4 triple junctions (!), but maybe it would square the system? I hope that's not required because the idea of doing these calculations on networks of complex faults with triple junctions and such would be amazing (e.g., rupture in quasi-fractal fault systems).

[mallickrishg]

Also, one random and unprincipled idea: What happens if we modify the model a little bit? Instead of just the three faults, what if we put the three faults inside a triangle that bounds the three interior elements? This would create 4 triple junctions (!), but maybe it would square the system? I hope that's not required because the idea of doing these calculations on networks of complex faults with triple junctions and such would be amazing (e.g., rupture in quasi-fractal fault systems).

Love the idea, but possibly too complicated for now. But eventually that is the sort of problem we will end up with if and when we tackle real branched/nested fault systems. For now, there might be a justifiable fix, but we will have to talk it out later today.

The conservation of slip at a point remains as a triple junction problem but the smoothness conditions can go back to same old 2 elements at a time equations - this makes the system square!

# triple junction equation
matrix_system[4,:] = matrix_slip[4,:] - matrix_slip[6,:] - matrix_slip[12,:]# x component
matrix_system[5,:] = matrix_slip[5,:] - matrix_slip[7,:] - matrix_slip[13,:]# y component
# regular smoothness
matrix_system[6,:] = matrix_slipgradient[4,:] - matrix_slipgradient[6,:]
matrix_system[7,:] = matrix_slipgradient[5,:] - matrix_slipgradient[7,:]
matrix_system[12,:] = matrix_slipgradient[4,:] - matrix_slipgradient[12,:]
matrix_system[13,:] = matrix_slipgradient[5,:] - matrix_slipgradient[13,:]

This is a different from

# smoothness at a triple junction
matrix_system[6,:] = matrix_slipgradient[4,:] - matrix_slipgradient[6,:] - matrix_slipgradient[12,:]# x component
matrix_system[7,:] = matrix_slipgradient[5,:] - matrix_slipgradient[7,:] - matrix_slipgradient[13,:]# y component

Both these equations seem to produce singularity free stress distributions. However, they are not the same! I don't even know which is right.

Triple junction smoothing

Back to regular smoothing (2 at a time)

[mallickrishg]

@brendanjmeade I have a new notebook in the exploratory/ folder that has the equations written out clearly for the triple junction implementation for both the methods we discussed today - https://github.com/brendanjmeade/bemcs/blob/main/exploratory/triplejunctions_boundaryconditions.ipynb

Let me know if there's something I need to clarify. Long story short - there continues to be a weak singularity at the triple junction. One of the methods is better than the other but they are not perfect.

[mallickrishg]

@brendanjmeade - seems like what we were missing was to explicitly specify boundary conditions on one of the nodes at the triple junction. However, it is unclear which of the three nodes we should specify BCs on. But the problem seems to have its own preference - which is baffling me! https://github.com/brendanjmeade/bemcs/blob/main/exploratory/triplejunctions_boundaryconditions.ipynb

[brendanjmeade]

@mallickrishg This notebook is so good! It's so clear and we can reason about it. I've read the part entitled "Recent fix" and it's an interesting way and certainly mystifying way to do this. When you write,

p3(n1)gradx = 0 or p2(n1)gradx = 0 (extra boundary condition)
p3(n1)grady = 0 or p2(n1)grady = 0 (extra boundary condition)

Exactly where are these being applied? On elements 3 or 2 but at which nodes? Intriguing!

[mallickrishg]

Thanks Brendan, I'm glad it makes sense. what I meant by those equations is that you can set slip gradient to 0 for any 1 out of the 3 nodes at the junction. But it seems to only lead to nonsingular stress when applied to either patch 2 or 3, but not to patch 1

• patch 1 node 3 slip gradient = 0 OR
• patch 2 node 3 slip gradient = 0

It is implemented in code as shown below,

# testing additional constraints
# matrix_system[12,:] = matrix_slipgradient[6,:]
# matrix_system[13,:] = matrix_slipgradient[7,:]
matrix_system[12,:] = matrix_slipgradient[12,:]
matrix_system[13,:] = matrix_slipgradient[13,:]

Regarding the actual equation number on the lhs (for matrix_system) I will come up with something more intuitive like - the equations will be stacked as all the equations at central nodes, then open nodes, then 2 overlapping elements and then triple junctions.

[brendanjmeade]

But it seems to only lead to nonsingular stress when applied to either patch 2 or 3, but not to patch 1

• patch 1 node 3 slip gradient = 0 OR
• patch 2 node 3 slip gradient = 0

That is something that I totally don't understand. I can see that adding and additional BC helps but I have to work to understand why it might be the case that there something special about why it only works for certainly patches. That feels like it has to be wrong because we can arbitrarily reorder the elements so that only other explanation would appear to be that there's some effect that depends on the geometry.

I'm reading through this now!

[mallickrishg]

I agree with your comment, it is so arbitrary it does not feel right. I think we should just stick with the normal formulation which is, 2 equations for continuity of slip at a triple junction and 4 equations for smoothness (2 nodes at a time).

I have now been able to come up with an automatic node labelling and linear operator construction algorithm. Could you look over it when you have time and tell me if the idea makes sense? It seems to work correctly - https://github.com/brendanjmeade/bemcs/blob/main/examples/examples_3_simpleintersectingfaults.ipynb

[brendanjmeade]

@mallickrishg

I've been working to improve and organize the basic triple junction work and have done so in:

https://github.com/brendanjmeade/bemcs/blob/main/exploratory/triple_junctions_boundary_conditions_bjm.ipynb

Here I have slightly modified your notebook with a core goal: Explore the different ways. I have started labeling these "cases" so that we can refer to them more easily. To enable this, I propose describing the system of equations in two parts: 1) the 12 rows that are the simple boundary conditions that we know and 2) 6 rows that we're unsure of. Each of the "cases" is a different way to implement these 6 rows. I think this is a clean way to document and explore these experiments.

I think we should continue exploring here, and I have two suggested todos:

• [ ] @brendanjmeade : Push some of the repetitive code into notebook local functions. These are not for permanence and so don't belong in bemcs.py but are intended to reduce repetitive code so that there are fewer chances for bugs to creep in.
• [ ] @mallickrishg : Rewrite matrix_system assembly. Put the 12 rows with the equations and boundary conditions we're confident in at the top of the matrix_system matrix. This will keep what we know well at the top of the matrix and have the six equations we're not sure of at the bottom.

The main goal of both of these changes is to get this code organized enough so that we can isolate exactly what makes sense and does not for each case. My current opinion is that we haven't got this right yet.

[brendanjmeade]

I agree with your comment, it is so arbitrary it does not feel right. I think we should just stick with the normal formulation which is, 2 equations for continuity of slip at a triple junction and 4 equations for smoothness (2 nodes at a time).

See previous comment. I don't think we've resolved the triple junction problem yet!

I have now been able to come up with an automatic node labelling and linear operator construction algorithm. Could you look over it when you have time and tell me if the idea makes sense? It seems to work correctly - main/examples/examples_3_simpleintersectingfaults.ipynb

I'm really excited about this but I want us to get the triple junction problem solved first. A lot hangs in the balance with that.

[brendanjmeade]

@mallickrishg Ok and it's late and I'm a bit tired but I have one more idea for the triple junction that requires relaxing one of our assumptions.

What if there aren't 12 boundary conditions that we think we know? What if instead we have this:

• For each of the three elements, specify BCs only at the open endpoint. That's six equations.
• At the triple junction:
  • Slip continuity (1-2, 1-3, 2-3) for x and y each. That's six equations.
  • Slip gradient continuity (1-2, 1-3, 2-3) for x and y each. That's six equations.

In other words, slip at the center of each element is not specified. There's no explicit notion of a triple junction in this case so I don't think slip will be properly conserved but I wanted to share that idea because it's a little bit different and is potentially useful for that reason.

[mallickrishg]

I love the new notebook that you made. It clearly describes the various options. In the new version of linear operator construction I am actually setting it up the way you suggest - all the known equations that we agree on at the top, and the triple junction equations at the bottom! Lots to try over the next few days

[brendanjmeade]

@mallickrishg

Maybe Case C is really close?!. Currently, the slip gradient continuity equations might be considered incomplete in that there are only two of them. One for the $x$ gradient of $u_x$ and one for the $y$ gradient of $u_y$. What if we just add two more? One for the $y$ gradient of $u_x$ and one for the $x$ gradient of $u_y$? That would give us 18 equations! I've written this out below.

Case C (Revised)

• Continuity of slip and slip gradient at a triple junction

bc#	expression	meaning
13	$p^x_1(n_3) - p^x_2(n_1) - p^x_3(n_1) = 0$	elements 1, 2, 3: continuity of $x$ slip at a triple junction
14	$p^y_1(n_3) - p^y_2(n_1) - p^y_3(n_1) = 0$	elements 1, 2, 3: continuity of $y$ slip at a triple junction
15	$\partial p^x_1(n_3) / \partial x - \partial p^x_2(n_1) / \partial x - \partial p^x_3(n_1) / \partial x = 0$	elements 1, 2, 3: continuity of $x$ slip gradient of $u_x$ at a triple junction
16	$\partial p^y_1(n_3) / \partial y - \partial p^y_2(n_1) / \partial y - \partial p^y_3(n_1) / \partial y = 0$	elements 1, 2, 3: continuity of $y$ slip gradient of $u_y$ at a triple junction
17	$\partial p^x_1(n_3) / \partial y - \partial p^x_2(n_1) / \partial y - \partial p^x_3(n_1) / \partial y = 0$	elements 1, 2, 3: continuity of $y$ slip gradient of $u_x$ at a triple junction
18	$\partial p^y_1(n_3) / \partial x - \partial p^y_2(n_1) / \partial x - \partial p^y_3(n_1) / \partial x = 0$	elements 1, 2, 3: continuity of $x$ slip gradient of $u_y$ at a triple junction

[mallickrishg]

I think this is a really interesting idea. The way I've been thinking about slip gradients has been to project the gradient tensor onto the fault shear direction to produce a gradient vector https://github.com/brendanjmeade/bemcs/blob/main/exploratory/derivation_slipgradient_smoothness.ipynb I need to think about it instead we should directly just work in tensor space and apply the smoothness conditions there, like you suggest?

[brendanjmeade]

Thanks for sharing the new notebook. There's good stuff here. I don't quite see where equation four comes from?

If we have a slip vector in some local ($s$, $n$) coordinate frame, we either need to rotate it or project it to some other frame that can be made common for all elements. The global frame ($x$, $y$) is an obvious candidate. Is equation 4 the projection via a dot product, as you do below for the slip gradient case? I like the way you wrote out and explained the slip gradient case.

After we get continuous smoothly varying slip at the TJ, I'll be interested to see if stresses are indeed continuous or just tractions.

[mallickrishg]

Regarding equation 4, i converted the slip vector from a local shear, normal coordinate system to an $x,y$ coordinate system. The conversion for any slip vector $\underline{s}$ from a $s,n$ coordinate system to $x,y$ coordinate system is below:

• $\underline{s_s}= |s_s| \hat{d}$ where $\hat{d}$ is the unit vector for shear in $x,y$
• $\underline{s_n}= |s_n| \hat{n}$ where $\hat{n}$ is the unit normal vector in $x,y$

[brendanjmeade]

I'm still not sure if I follow this but I'll sit with it for a while and see. I'd usually write this as $\mathbf{s}_G = \mathbf{R} \mathbf{s}_L$, where $\mathbf{R}$ is some rotation matrix unique for each element to takes us from local ($L$) to global ($G$) coordinate systems.

[mallickrishg]

I think what we have should be identical. The rotation matrix for me is something like $R = \begin{bmatrix} d_x & n_x \ d_z & n_z \end{bmatrix}$

[mallickrishg]

@brendanjmeade I am going through your triple_junctions_boundary_conditions_bjm.ipynb and I thought I can answer some questions you had in your notes

Case A

• Does slip actually vary smoothly from element to element at the triple junction? I'm asking this because the displacement as a function x figure looks non-smooth at the triple junction. Also equation 15-18 cover the 1-2 and 1-3 smoothness case but what about the 2-3 case? We could write the 2-3 case but I think you mentioned that this would be redundant?

This is a really interesting thing - the triple junction slip condition does not make slip continuous if you only look at only 2 elements because it partitions slip from element 1 to elements 2 + 3, so slip is still conserved. Slip and as a consequence displacement ends up being non-smooth and as a result the stresses can be discontinuous. However, I have been feeling comfortable with this solution because even though stresses are discontinuous, they are not singular, and when you project the stress tensor to the appropriate traction vector, it looks smooth. I did not do this very carefully because it requires a bit more work to define the unit normals correctly in space to project the stress tensor and visualize spatially varying traction vectors along multiple fault planes. Perhaps this is something I need to do at some point.

Case B

• Continuity of slip at a triple junction using element by element pairing. This is essentially what we do for the two element case and there is no direct notion of a triple junction in the math here.

You are absolutely correct. These stresses are singular! It's basically what we would get if we had linear basis functions - spatially continuous but not differentiable slip.

Case C

• Continuity of slip at a triple junction
• Smoothness of slip between any 2 elements that share a common node

Once again you are right - it is a bizarre boundary condition indeed! In the example in your notebook, the reason why those stresses look so good (continuous) is because we have set the slip gradient for the triple junction node of element 3 to 0 and only plotted stresses along a transect that traverses element 1 to element 2. Now if we changed the geometry such that element 2 was vertical and element 3 was horizontal, you will see stresses appear discontinuous. But if you now set the slip gradient for element 2 to 0, again the stresses will look continuous along this transect. The condition number for these 2 interchangeable cases is almost identical. However, now if we apply the slip gradient of element 1 = 0, everything falls apart. Large condition number and stresses are singular. But that does not make sense to me!

I wonder if the reason for this issue is that the slip gradient condition is $\nabla s_1 - \nabla s_2 - \nabla s_3 = 0$ and the extra BC is $\nabla s_3 = 0$ or $\nabla s_2 = 0$ or $\nabla s_1 = 0$ which is equivalent to: $\nabla s_1 - \nabla s_2 = 0$ and $\nabla s_3 = 0$ with interchangeable $s_2$ and $s_3$ However for $s_1$, the equations might need to be $\nabla s_2 - \nabla s_3 = 0$ and $\nabla s_1 = 0$? If this is the case, it means all of these choices are 'correct' in some sense, but return different results of partitioning of slip between the 3 elements.

Case D

My hope is that we can completely avoid this if we switch to working on the tensor equations for the Jacobian of slip. I can explore this tomorrow.

[brendanjmeade]

@mallickrishg Thanks for the great responses. With regard to Case C I think you are right that adding a gradient constraint on $s_1$ essentially introduces a type of linear dependence that leads to the degeneracy. I'm hoping that the full vector version fixes that!

[mallickrishg]

I spent some time thinking about the triple junction problem and how we have so many options to solve the governing equations. I realized that some of the options don't make as much sense if we were to generalize to M-elements at a junction! M fault elements at a junction requires 6M equations. 2M of these equations are provided at central nodes and 2M at open nodes, leaving us with 2M equations that we need to close the system.

• CASE D: Conserving the gradient of slip in a tensor sense completely fails when the number of elements at a junction exceeds 3. It results in too few equations (it can only provide 4 gradient equations) and the system is underdetermined! This approach also fails when we consider only 2 fault elements - then we have too many equations. It seems to only make sense for a triple junction (even in that case, it might work only in theory - have not tested this).
• CASE B: Since this only has continuity of slip for 2 elements at a time, we end up with 2M equations to close the system but only 2(M-1) of them are unique AND the stresses are singular
• CASE C: We apply 2 equations for conservation of slip at the junction, 2 equations for continuity of slip gradients across 2 arbitrarily chosen elements and 2(M-2) equations setting the slip gradient to 0 for the remaining M-2 elements. This case will presumably continue to produce non-singular stresses but it seems really really arbitrary now for the general M-element case.
• CASE A: possibly the most defensibly case where we provide 2 equations for conservation of slip at the junction, and 2(M-1) equations equating the slip gradient for M-1 pairs of nodes. This closes the system and does not have any arbitrary choices for which elements which conditions are applied to. Currently the condition I am strongly leaning on, and what I propose as our panacea.

Once the rewrite and testing is done, I intend to test this M-element at a junction with case A to convince both of us that we have solved the problem.

[brendanjmeade]

I hadn't even considered dealing with junctions with $M>3$. I think that if we have a well-posed solution to the 3-junction case (hopefully the matching derivatives approach will do this!) then we might consider only having triple junctions and just putting two or more really close to each other if we wanted to emulate an $M>3$ case. That would also avoid us introducing a structure that is inherently kinematically unstable (at least in the plate tectonics sense). I'll think on this more!

[mallickrishg]

I don't intend for us to deal with junctions of elements $M>3$. As you noted, they are not kinematically stable and may not be really meaningful scientifically (at least in the fracture and earthquake mechanics communities - IMO).

However, I was using the general $M$ element case to argue that most of our methods of dealing with triple junctions are rendered unusable except for case A i.e., 2 conservation of slip equations and 2(M-1) slip gradient equations.

[brendanjmeade]

I hadn't appreciated that perspective!

[mallickrishg]

@brendanjmeade good news! the jacobian of slip actually works quite nicely. I have added it as a section to the existing triple junction notebook - https://github.com/brendanjmeade/bemcs/blob/main/refactor/example_triple_junctions.ipynb

It gives discontinuous stresses, but they are NOT singular. In fact, for many combinations of slip and fault geometries I found that the displacements and stresses at the triple junction are quite close (not identical though) for case A (my preferred case that can be generalized for N elements) and case D (jacobian of slip). Attached is one such test case.

Just to remind you: case A: conservation of slip (2 equations), directional gradient of slip (2 elements at a time) are continuous (4 equations) case C: conservation of slip (2 equations), directional gradient of slip for 2 selected elements is continuous, while the third element's gradient is set to 0 case D: conservation of slip (2 equations), Jacobian of slip is conserved (4 equations)

[brendanjmeade]

That's a really interesting experiment and thanks for pushing on this. Notably, all the solutions look pretty similar, which is something that I hadn't appreciated or anticipated. I'm still not sure what to make of this. Does it mean that the BCs essentially don't matter? Are the tractions continuous in all 3 cases? To help me better understand this, I'll start plotting the volumetric displacements and stresses tomorrow.

[brendanjmeade]

@mallickrishg

I'm starting work on the triple junction notebook a bit, and I had a question about how slip is specified on the three elements. I see this written as:

slip_vector_x = np.array([1.0, 0.8, -0.0])
slip_vector_y = np.array([0.5, 0.5, -0.5])

I see they x and y labels that strongly suggest that these are in a global Cartesian reference frame. I'm writing to check on that. The reason I'm asking is that I'm looking to build some examples with specific slip distributions.

[mallickrishg]

@brendanjmeade slip is indeed specified in a global coordinate system through matrix_slip as the linear operator

[brendanjmeade]

@mallickrishg Got it...and then they get rotated in fault local slip components in get_matrices_slip_slip_gradient in bemcs.py?

[mallickrishg]

@brendanjmeade I just realized that when I compute tractions on the fault, they are not continuous for the triple junction. Then I went back and looked at 2 non-planar elements and even there the tractions have small discontinuities. I think the continuity of slip and smoothness equations only guarantee that singularities will not exist but I think discontinuties in stress & tractions don't seem to go away.

My guess is that it is related to a point you had brought up a while back about C0 and C1 meshes. Because our mesh is not smooth, we might never be able to make tractions smooth.

Does this affect strain energy estimates? Perhaps, but the discontinuous nature of stress does not mean any of these quantities are unbounded. So, I think we are still fine.

[mallickrishg]

@mallickrishg Got it...and then they get rotated in fault local slip components in get_matrices_slip_slip_gradient in bemcs.py?

get_matrices_slip_slip_gradient rotates the 3qn coefficients, which are in local (shear, normal) coordinates to a global (x,y) coordinate system. So both LHS and RHS end up in (x,y).

[brendanjmeade]

Then I went back and looked at 2 non-planar elements and even there the tractions have small discontinuities. I think the continuity of slip and smoothness equations only guarantee that singularities will not exist but I think discontinuties in stress & tractions don't seem to go away.

I hadn't noticed this before and am surprised because I would have thought this would be very easy to see/identify. Do you have a notebook the shows this? This information suggests that something is wrong here. I say this because there's nothing discontinuous (normal vector is the same for all elements) and slip is supposed to be continuous everywhere here except at the endpoints. Thinking quickly, I can see a few things that we might be doing wrong.

• The quadratic Greens' functions are incorrect.
• There is some mixup/ambiguity related to applying slip continuity at element end points and then translating that to what the coefficient weights should be at the three quadratic element nodes, which do not include the endpoints.

Right now I think this is a problem. What do you think?

[mallickrishg]

I say this because there's nothing discontinuous (normal vector is the same for all elements) and slip is supposed to be continuous everywhere here except at the endpoints.

It's only an issue for non planar faults - the normal vector for the 2 elements are not the same at the overlapping node. However slip and the slip gradient are continuous, so it's a mystery to me as to why the tractions are discontinuous. Here is the notebook where I tried it out - https://github.com/brendanjmeade/bemcs/blob/main/notebooks/example_non_planar_fault.ipynb

[brendanjmeade]

I see and this makes sense. So everything is still perfectly continuous for planar faults right? I think I misread your prior email!

[mallickrishg]

I need to check that planar example thoroughly too. It turns out the last time I made this plot for the non-planar case, the slip vector values were just right that there wasn't a discontinuity in tractions and so I was fooled.

[mallickrishg]

just checked planar faults with non-uniform lengths have no problem with continuity of stress: https://github.com/brendanjmeade/bemcs/blob/main/notebooks/example_non_uniform_planar_patches.ipynb