Position-based simplex quadrature?

[zasexton]

Context

Discussion #587 and Issue #588 identified that the Niederer CEP benchmark supplied Quadrature_modifier_TRI3 and Quadrature_modifier_TET4 values even though its volume mesh contains only HEX8 elements and its boundary meshes contain only QUAD4 elements. Removing those settings corrects the benchmark, but it leaves a broader design question: should the position-based quadrature controls remain part of the public input interface, and if so, how should they be represented and documented?

The controls were introduced in PR #56 and Issue #59 as non-Gaussian quadrature options intended to reduce conduction-velocity sensitivity to mesh resolution. These rules can be calibrated to reduce conduction-velocity variation for a specified problem and mesh-size interval.

What the current controls do

Both controls move equal-weight quadrature points along centroid-to-vertex lines while preserving a symmetric rule:

• Quadrature_modifier_TRI3 is attached to an Add_face, accepts a value s in [1/3, 1], and uses the three barycentric permutations of (s, t, t), where t = (1-s)/2. Its Gaussian default is s = 2/3.
• Quadrature_modifier_TET4 is attached to an Add_mesh, accepts a value s in [1/4, 1], and uses the four barycentric permutations of (s, t, t, t), where t = (1-s)/3. Its Gaussian default is s = (5 + 3*sqrt(5))/20, approximately 0.5854102.

The weights remain fixed. Every accepted value preserves exact integration of constants and linear polynomials. Within each accepted family, however, only the Gaussian default is exact for quadratic polynomials. A nondefault position therefore changes the rule from degree-two to degree-one polynomial exactness.

This also means that selecting only an "exactness level" would not fully represent the existing feature. All nondefault positions have the same degree-one exactness, but different positions can produce different numerical dispersion and conduction velocities. An exactness-only interface could select a conventional reduced rule, but it may not reproduce a position calibrated in a published study from @oguzziya2.

There are additional scope questions in the current interface:

• The TRI3 setting is a face-level control, whereas the TET4 setting is a mesh-level control.
• The TET4 setting is not scoped specifically to CEP, even though its motivating use case was cardiac electrophysiology.
• A setting is accepted even when the corresponding element type is absent, which allowed the Niederer benchmark inputs to silently advertise behavior they did not exercise.
• The input names, coordinate convention, fixed weights, defaults, ranges, and change in polynomial exactness are not explained in the User Guide.

What the published examples establish

Woodworth, Cansız, and Kaliske (2022) describe modified quadrature as placing points at arbitrary isoparametric offsets. Their numerical study examines dependence on spatial discretization and conductivity and demonstrates improved conduction-velocity accuracy on selected coarse meshes.

Tikenoğulları et al. (2023) use the method for a growth-and-remodeling problem in which the mesh deforms. They identify the quadrature locations giving the flattest conduction-velocity response over the element-size ranges expected between their baseline and grown models. Their calibration is therefore tied to a particular mesh-size interval and electrophysiology setup rather than establishing invariance for arbitrary resolutions, mesh families, conductivities, time steps, or propagation directions.

Note: In reviewing the paper, it calibrates line and tetrahedral rules using its own quadrature-location convention, while the current svMultiPhysics inputs expose TRI3-face and TET4-mesh parameters named only as "modifiers." Without an explicit definition of the coordinate mapping, users should not assume that a numerical value can be copied directly from a paper into an XML input.

Options

1. Remove or deprecate the position controls

Use only conventional quadrature rules. This gives the clearest accuracy contract but makes it harder to reproduce prior studies that intentionally calibrated quadrature-point locations.

2. Retain the existing parameters and improve their documentation

This preserves input compatibility, but the generic modifier name would continue to hide what is being changed and why a nondefault value is risky.

3. Expose an explicit quadrature-rule choice

Make standard Gaussian quadrature the default and require users to explicitly select a position-based custom rule before supplying s. Conceptually, the interface would distinguish among:

• standard Gaussian quadrature;
• central or nodal rules, if these are intentionally supported; and
• a position-based symmetric rule with a documented vertex-associated barycentric coordinate.

This option preserves research reproducibility without presenting a calibrated position as a general-purpose accuracy setting. Polynomial exactness should be documented and reported as a property of the selected rule, not used as the only selector.

Design Questions:

• Are there current svMultiPhysics models, beyond the cited studies, that rely on nondefault TRI3 or TET4 values?
• Should TRI3-face and TET4-mesh controls be considered separately because they have different scopes and use cases?
• Is backward compatibility with the exact scalar controls necessary, or would a deprecation period be allowable?
• Would a documented calibration example and validation dataset make this feature maintainable, or should it be removed until such evidence exists?

[zasexton]

Additionally, how should position-based quadrature calibration behave for domains containing different element topologies? The current modifier is specific to TET4 volume elements, while other supported topologies retain their standard quadrature rules. This can make numerical conduction velocity depend on element topology and introduce inconsistent propagation error across topology transitions, even when nominal mesh resolution and material parameters are unchanged. The design should either define topology-specific calibration and verification for every supported element type, explicitly restrict the feature to homogeneous TET4 meshes, or avoid presenting it as a general mesh-independence mechanism. Mixed-topology input should also be validated or partitioned explicitly into topology-specific element blocks.

[michelebucelli]

Thank you @zasexton for looking into this!

I am not familiar with these modified quadrature rules, so I have no strong opinion on their usefulness. Based on the information you and @javijv4 retrieved from those papers, I agree with your take that the modified formulas should not be the default behavior, but rather that standard Gaussian quadrature should.

Concerning how to implement and expose this feature to the user, of the options you propose, I think 3 (expose an explicit quadrature rule choice) is the best.

I do not know what your exact plans are for refactoring the quadrature formulas, but in my mind the minimal interface for a quadrature formula (and thus the most abstract class) is a list of weight-point pairs (this is what deal.II does, for example). From the abstract interface you could derive a concrete class for Gaussian quadrature, and separately a concrete class for these modifed quadrature formulas. The former would take in input (from the constructor, or from parameters) the desired degree of exactness, while the latter would take in input the shift value (and maybe the desired degree of exactness, not sure if it makes sense in this context). The rules would/should be similarly distinct on the user-facing side (i.e. in the XML file). I think this would minimize any risk of misunderstanding.

[javijv4]

@zasexton, I also think option 3 is the best, as it keeps the default to the standard Gaussian rules and makes it clear to the user whether they want to use the modified ones.

As far as I am aware, all papers from @MatteoSalvador used these rules and @elenasmartinez last publication:
* Martinez2026
* Salvador2026
* Salvador2024a
* Salvador2024b

I do think it is important to do the topology check, such that the TET4 modifier can only be used in linear tetrahedral meshes. Same for TRI3.

I also think it would be ideal if there were an example or validation study. Maybe using the Niederer benchmark with a tet mesh would be enough? But I also wonder if @MatteoSalvador has any thoughts on this.

[zasexton]

I will proceed to add notes for this option 3 route to the currently open quadrature refactor issue #580 . Ultimately, I think it might be nice to sit down with the CEP people in the lab to understand more about how and when this integration rule is sufficient. The paper from Oguz seems to optimize the modified quadrature points to minimize CV mesh sensitivities. I don't know if these modified non-gaussian quadratures also guarantee convergence accuracy since I don't believe the original paper provides numerical convergence evidence for the modified values. I'll plan to have this as a separate explicit rule.

[ktbolt]

@zasexton Option 3 looks like the best choice to me.

[MatteoSalvador]

@zasexton apologies for the late response, last week has been really busy for me.

To give some additional background, these modified quadrature rules have been thoroughly validated and tested directly on patient-specific HLHS ventricular geometries by Oguz and myself (including convergence studies). Since 3+ years ago running EP on a very fine TET meshes (i.e. mesh size approximately equal to 0.35 mm or less) in (back then) svFSIplus would have been really inefficient from the computational perspective, this approach was the only way to obtain reasonable/realistic conduction velocities on coarser meshes (between 0.5 mm and 1.5 mm in terms of average mesh size according to the case at hand).

I agree that the CI tests involving Niederer benchmark with HEX elements + these quadrature modifiers are misleading and should instead run on a tetrahedral version of the benchmark to properly cover these code paths. This was a mistake on my side.

In terms of design options moving forward, I agree with @ktbolt that option 3 seems to be the best choice.

I also want to highlight that during my postdoc I thoroughly validated and tested patient-specific EP pipelines involving Monodomain equation coupled with TTP. For this reason, I do not exclude that there might be issues arising with other ionic models (e.g. see #540) or other code paths that I did not consider in my research papers (e.g. HEX vs. TET meshes).

Of course, feel free to ping me if you observe any unexpected behavior in the pipelines described in my publications (@javijv4, @elenasmartinez and @michelebucelli). The library evolved a lot in the last 2 years and I lost track of the changes that happened after my postdoc.

[zasexton]

I have proceed to set up a tetrahedral Niederer benchmark examining Gaussian spatial and temporal controls with two target-mesh tetrahedralizations. I examine the modified TET4 s=0.34 rule, and a sweep of the TET4 quadrature position at Δx=0.1 mm and Δt=0.005 ms. Some preliminary findings are:

• The historical non-Gaussian TET4 modifier s=0.34 value was less accurate than Gaussian integration at the target resolution. Relative to the high-resolution Gaussian comparison, activation-field RMSE increased from approximately 0.373 to 1.140 ms, far-corner error increased from 1.041 to 2.419 ms, fiber-direction apparent-CV error increased from 0.441% to 1.876%, and transverse-direction apparent-CV error increased from 0.629% to 7.597%.
• The best modifier depends on the output being matched. The sweep favors approximately s=0.60 for transverse apparent CV, s≈0.63 for fiber apparent CV, and a larger value near s≈0.68 for far-corner activation.
• The best all-around non-Gaussian modifier candidate, s=0.65, reduced activation RMSE to approximately 0.159 ms and fiber apparent-CV error to 0.178%, but increased transverse apparent-CV error to 2.086%.

It appears that activation-field error, fiber apparent CV, and transverse apparent CV favor different quadrature positions.

Additionally, the far-corner activation was non-monotonic over the spatial refinement, and the fiber apparent-CV produced a negative observed convergence order. The 0.05 mm / 0.001 ms result was 41.4255 ms at P8, compared with the 42.53–43.18 ms refined results reported by the Niederer Benchmark study. Although the published results here are an aggregate of implementations and not necessarily a formal acceptance window.

I would still need to run additional high-resolution Gaussian quadrature simulations to confirm this trend since the current high-res reference Gaussian simulation (Δx=0.05 mm and Δt=0.001 ms) remains unqualified. The the Gaussian results did not converge consistently as the mesh was refined—the activation times changed direction, and the fiber conduction-velocity differences grew instead of shrinking. Consequently, the finest 0.05 mm / 0.001 ms result cannot show that it is sufficiently close to mesh- and time-step-independent solution sensitivity.

Thank you @MatteoSalvador for the additional context! Generally, I am just wanting to determine how transferable non-Gaussian rules might be treated to inform new users that might be unfamiliar with such Gaussian modification parameters. Specifically in what parameter ranges we can be comfortable documenting mesh-resolution-independence settings for the code as it is refactored.

I'll continue to run a couple more high-resolution Gaussian models to see if we can reach a qualified simulation to compare against.