🐛[BUG]: total surface area can be wrong with stl tessellation samples

[Karl-JT]

Version

1.3.0

On which installation method(s) does this occur?

Docker

Describe the issue

When using an stl file (sphere for example) instead of primitive geometry, the total area of all sampled points are different from the actual total area. This leads to errors on computing drag or lift when summing up the pressure on the surface. The reason is in tessellation geometry sampling. The area of samples are calculated based on each triangle, but some triangle will have zero samples.

Minimum reproducible example

@modulus.sym.main(config_path="conf", config_name="config")
def run(cfg: ModulusConfig) -> None:
    point_path = to_absolute_path("./stl_files")
    sphere_geo = Tessellation.from_stl(
        point_path + "/sphere.stl", airtight=True
    )

    ns = NavierStokes(nu=1.0, dim=3, time=False)

    # determine inputs outputs of the network
    input_keys = [Key("x"), Key("y"), Key("z")]
    output_keys = [Key("u"), Key("v"), Key("w"), Key("p")]

    flow_net = FullyConnectedArch(
        input_keys=input_keys,
        output_keys=output_keys,
        adaptive_activations=cfg.custom.adaptive_activations,
    )

    nodes = (
        ns.make_nodes()
        + [flow_net.make_node(name="flow_network")]
    )

    domain = Domain()

    interior = PointwiseInteriorConstraint(
        nodes=nodes,
        geometry=sphere_geo,
        outvar={"continuity": 0, "momentum_x": 0, "momentum_y": 0, "momentum_z": 0},
        batch_size=1000,
        lambda_weighting={
            "continuity": Symbol("sdf"),
            "momentum_x": Symbol("sdf"),
            "momentum_y": Symbol("sdf"),
            "momentum_z": Symbol("sdf"),
        },
    )
    domain.add_constraint(interior, "interior")

    force_p = PointwiseMonitor(
        sphere_geo.sample_boundary(100),
        output_names=["p"],
        metrics={
            "total_area": lambda var: torch.sum(var["area"]),
        },
        nodes=nodes,
    )
    domain.add_monitor(force_p)

    # make solver
    slv = Solver(cfg, domain)

    # start solver
    slv.solve()


if __name__ == "__main__":
    run()

Relevant log output

RuntimeWarning: divide by zero encountered in divide
  np.full(x.shape, triangle_areas[index] / x.shape[0])

The monitored total surface is 9866.225, but the actual surface of the sphere is 31415.93.

Environment details

No response

Other/Misc.

No response

[Karl-JT]

One potential solution could be change the following

invar["area"].append(
    np.full(x.shape, triangle_areas[index] / x.shape[0])
)

to

invar["area"].append(
    np.full(x.shape, np.sum(triangle_areas) / nr_points)
)

or to enforce at least one sample for each triangle.