[BUG] LoopProjectfileProcessor hardcodes fault_stratigraphy=None, producing models with no fault offset

[tamaravasey]

Describe your issue

Summary

LoopStructural.modelling.input.project_file.LoopProjectfileProcessor.__init__
passes fault_stratigraphy=None to ProcessInputData.__init__. When
GeologicalModel.from_processor(processor) then builds each foliation, it
calls create_and_add_foliation(s, faults=None) and the stratigraphy feature
ends up with zero fault references on its builder. The implicit foliation
field is built smoothly across every fault and extracted isosurfaces are
continuous meshes that show no fault offset — even though the fault features
themselves are interpolated and exist in the model.

The data needed to construct the proper mapping is already in the .loop3d
(eventRelationships carries the unit-fault crossings that map2loop's
Stage 4 topology produces). It is simply never read.

The result: the same map2loop output produces two different LoopStructural
models depending on which processor consumes it. Map2LoopProcessor
builds fault_stratigraphy correctly from its fault_strat DataFrame;
LoopProjectfileProcessor discards the relationships. The .loop3d is
supposed to be a faithful re-encoding, and today it is lossy for the
unit-fault mapping.

Versions affected

	Version	fault_stratigraphy=None hardcoded
LoopStructural released	1.6.27	✓
LoopStructural released	1.6.28	✓
LoopStructural master HEAD	(current)	✓

Confirmed on Python 3.12, conda-forge / loop3d channels. LoopProjectFile
0.2.2, map2loop 3.3.1.

Reproduction

from LoopProjectFile import ProjectFile
from LoopStructural.modelling.input.project_file import LoopProjectfileProcessor
from LoopStructural import GeologicalModel

# Any .loop3d with stratigraphy + faults works. E.g. the bundled map2loop
# Hamersley demo (Project(use_australian_state_data="WA", ...).run_all()).
pf = ProjectFile("hamersley.loop3d")
proc = LoopProjectfileProcessor(pf)

assert proc.fault_stratigraphy is None  # the bug — should be a dict

model = GeologicalModel.from_processor(proc)
model.update()

sg = model['sg']
print(len(sg.faults))  # 0 — should be N

The foliation has zero fault references even though the model knows
about every fault.

Empirical evidence

Sampling the foliation scalar field on a horizontal slice of the bundled
Hamersley demo, both with the existing None and with
_fault_stratigraphy = {sg: [all_fault_names]} set on the processor:

Configuration                           | mean |Δfield| | nonzero diff cells
----------------------------------------+--------------+-------------------
fault_stratigraphy=None  (current)      |   reference  |   reference
fault_stratigraphy populated, 1× disp   |        0.00  |   0 / 15717
fault_stratigraphy populated, 50× disp  |     2248.91  |  15717 / 15717

The 1× test is identical to within float precision — not because the fix
is a no-op, but because map2loop emits a placeholder 100 m displacement
when the source has no throw measurements, and 100 m is invisible against
the Hamersley regional ~6 km vertical signal. The 50× test (5000 m,
geologically realistic) produces visibly chopped layers and a synclinal
warp — proving the wiring works end-to-end.

Asymmetry between processors

	Map2LoopProcessor	LoopProjectfileProcessor
Source	map2loop output dir	.loop3d (re-encoded map2loop output)
Builds fault_stratigraphy	✓ from fault_strat DataFrame	✗ hardcodes None
Stratigraphy gets fault refs	✓	✗
Visible fault offset in surfaces	✓	✗

Root cause

LoopStructural/modelling/input/project_file.py, in
LoopProjectfileProcessor.__init__:

        super().__init__(
            ...
            fault_properties=fault_properties,
            fault_edges=fault_edges,
            colours=colours,
            fault_stratigraphy=None,      # ← here
            intrusions=None,
            ...
        )

Secondary: fault_stratigraphy is exposed on ProcessInputData
as a property without a setter, so users who notice cannot set it
post-construction through the public API — they must write to
_fault_stratigraphy directly.

Proposed fix

Two-tier, with a small extra:

Tier 1 — read the relationships from the .loop3d.
The data is in pf.eventRelationships (or the equivalent eventLinks /
faultLog tables Stage 4 topology populates). Build the mapping the same
shape Map2LoopProcessor produces:

# Sketch — exact column names depend on the .loop3d schema
fault_stratigraphy = {sg: [] for sg in stratigraphic_supergroups}
for row in pf.eventRelationships.itertuples():
    fname = self._fault_id_to_name.get(row.eventId1)
    uname = self._strat_id_to_name.get(row.eventId2)
    if fname and uname:
        fault_stratigraphy[supergroup_for(uname)].append(fname)
fault_stratigraphy = {sg: sorted(set(f)) for sg, f in fault_stratigraphy.items() if f}

Tier 2 — fallback. When eventRelationships is empty/missing,
default to "every fault cuts every supergroup":

if not fault_stratigraphy:
    fault_stratigraphy = {sg: list(fault_names) for sg in stratigraphic_supergroups}

Strictly better than None, and matches the standard assumption in
LoopStructural tutorials when users don't supply per-unit relationships.

Tier 3. Add a setter to the fault_stratigraphy property on
ProcessInputData so consumers can set it post-construction.

Workaround until the fix lands

processor = LoopProjectfileProcessor(pf)
processor._fault_stratigraphy = {  # write to backing field; property has no setter
    sg: list(processor.fault_network.faults)
    for sg in processor.stratigraphic_column.keys()
    if sg != "faults"
}
model = GeologicalModel.from_processor(processor)
model.update()

Applies Tier 2 only (every fault cuts every supergroup). Produces correct
fault wiring; visible offset still depends on having realistic displacement
values in the source data — map2loop's default 100 m placeholder is too
small to resolve at regional scales.

Why this matters

Regional models built from .loop3d files have been silently shipping
without fault offset on the stratigraphy. Users seeing layered models
with floating fault planes (rather than chopped, offset stratigraphy)
likely interpreted it as a data-quality issue, not a wiring bug. The
reference Brockman Syncline model at
tectonique.net/models/brockman_syncline.html
clearly shows fault offset — making the bug visible by comparison whenever
a user runs the bundled demo and compares.

PR offer

Happy to submit a PR implementing Tier 1 + Tier 2 + the setter, with a
regression test that builds a small model from a .loop3d and asserts
len(model['sg'].faults) > 0.

Minimal reproducing code example

In the text above

Error message

Note above

[lachlangrose]

Hi @tamaravasey thanks for opening this issue, happy for you to make a pr and I can review it as I don't know when I'll have a chance to fix the issue.

FYI I'm working on a new version of loopstructural so the api will change significantly especially around the model format and design. It will probably become live in a month or so, so depending on your use cases it might be worthwhile waiting for that.