Ordering of fault elements needs to be updated

[mallickrishg]

@brendanjmeade - i noticed there is an issue with how we deal with overlapping nodes for 2 fault elements. I don't have a fix off the top of my head, but I can think about it. Solving this problem is not crucial currently but it is something we need to tackle as we work with finite domains and closed polygons. Below I describe the issue and show 2 figures.

For 2 horizontal fault elements, bemcs.standardize_els_geometry handles the element ordering issue as follows:

if els.x2[i] < els.x1[i]:
   els.x2[i], els.x1[i] = els.x1[i], els.x2[i]
   els.y2[i], els.y1[i] = els.y1[i], els.y2[i]

so the unit vector for the 2 elements are forced to be in the same direction and the slip calculations work out as we expect.

However, if we were to break this organization of end points, and come up with an example of a vertical fault - I show below how the solution goes haywire when we changed the ordering of end points. I even use my automatic labelling scheme - and looks like there's something incorrect in the way we are dealing with the quadratic coefficients and fault end points.

Thoughts?

This is the "regular" case - displacements and slip are continuous and makes sense

This is the case of "flipped coordinates" for 1 element, and the linear operator is singular! This is a solution with pinv!

[mallickrishg]

The notebook containing this example is here: https://github.com/brendanjmeade/bemcs/blob/hrf-bem/examples/example_two_non_planar_fault_elements.ipynb

[brendanjmeade]

However, if we were to break this organization of end points, and come up with an example of a vertical fault - I show below how the solution goes haywire when we changed the ordering of end points. I even use my automatic labelling scheme - and looks like there's something incorrect in the way we are dealing with the quadratic coefficients and fault end points.

Thoughts?

I think the current node ordering approach work as long as faults aren't vertical. My hacky suggestion is that we should simply perturb vertical fault so that they are nearly vertical. There is a similar challenge in the block closure algorithm that we deal with by offsetting one end point by some small value (e.g., 1e-10). Would that be sufficient of or are there other cases too?

[mallickrishg]

I think what you suggest works. However, the problem that I foresee is on dealing with input files where this reordering is not possible i.e., consider the mht input file: if the right edge is meshed from bottom to top and topo is labelled from left to right, then the location where they meet can end up with unit normals in the "wrong" order, and we would have to reorder all the nodes for either right_edge or topo. Basically, we need to ensure the circulation within a polygon is in the same direction - which as we've discussed before is not an easy problem.

I'm trying to figure out how to deal with the quadratic coefficients for overlapping nodes agnostic to the circulation pattern of the elements. I'll keep you posted on if it works.

[brendanjmeade]

My suggestion is to avoid purely vertical segments at the input file stage. We can write a little checking function/script that looks for this and tells us which elements are invalid because they are vertical, and leave it to the user to fix the input file. That's probably a solution that is and efficient use of our time at this point. Thoughts?

[mallickrishg]

Unfortunately this is not an issue with vertical faults, the issue is in how we order the quadratic coefficients relative to end points. Even for acute angle fault intersections (see figure), when we reorder the $x,y$ coordinates, this problem persists.

[mallickrishg]

I'm pursuing this because I worry that when we work on more complex geometries, this book-keeping is going to be one of the more complicated parts of our task. See figure in previous post for the HRF problem.

But good news is that all we really need to do I think is to get coefficients for the second element (the one with 'inward' facing unit vector) to be reversed in order. I ran some simple tests, if i just flip the order of shear and tensile coefficients - i get the correct results! Now to figure out how to setup the linear operator to give us that!

[brendanjmeade]

Ahh I see and understand. Thanks for the examples! I looked back at Bem2d.jl where I tried to do a model with two boxes. There are two notebooks for this:

• https://github.com/brendanjmeade/Bem2d.jl/blob/master/examples/layeredboxfault.jl
• https://github.com/brendanjmeade/Bem2d.jl/blob/master/examples/layeredboxfault02.jl

I clearly have notes in these files about flipping normals...but I don't know any of it ever worked? I definitely agree that this is a very important problem to solve, and perhaps one of the most critical so that we can do material property variations. I honestly really want to work on this problem right now but I'm putting my bemcs effort into writing because I want the first paper out the door. Can't wait to get that done and am making good progress.

[mallickrishg]

Thanks for sharing and for working on the draft. At this stage, I think it's fine that we haven't solved everything. At least we know how to deal with mesh/ node issues if they come up. I realized the source of those near infinite condition numbers was related to the flipping of normals - we just need to order the mesh carefully.

[mallickrishg]

Closing this issue because we have a temporary fix for the problem - I have added flags/checks to deal with problematic nodes: https://github.com/brendanjmeade/bemcs/blob/main/bemcs/bemcs.py#L2260 Will need to figure out issues with circulation of mesh elements another day!