From e27baef92ddca7383746acc79aaebde1f3c61a9d Mon Sep 17 00:00:00 2001
From: Kaustubh Tangsali <71059996+ktangsali@users.noreply.github.com>
Date: Thu, 14 Nov 2024 22:42:22 -0800
Subject: [PATCH] [Bug] Fix area calculations for STLs (#205)

* fix the area calculations

* Increase tolerance
---
 CHANGELOG.md                             |  2 ++
 modulus/sym/geometry/tessellation.py     | 40 ++++++++----------------
 test/test_pdes/test_linear_elasticity.py | 12 +++++--
 3 files changed, 24 insertions(+), 30 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 89515b24..955dec65 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -22,6 +22,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Fixed
 
+- Fix the area calculated for STL meshes.
+
 ### Security
 
 ### Dependencies
diff --git a/modulus/sym/geometry/tessellation.py b/modulus/sym/geometry/tessellation.py
index ecd4252a..b87c69cb 100644
--- a/modulus/sym/geometry/tessellation.py
+++ b/modulus/sym/geometry/tessellation.py
@@ -95,16 +95,17 @@ def sample(
                     invar["normal_z"].append(
                         np.full(x.shape, mesh.normals[index, 2]) / normal_scale
                     )
-                    invar["area"].append(
-                        np.full(x.shape, triangle_areas[index] / x.shape[0])
-                    )
+
                 invar["x"] = np.concatenate(invar["x"], axis=0)
                 invar["y"] = np.concatenate(invar["y"], axis=0)
                 invar["z"] = np.concatenate(invar["z"], axis=0)
                 invar["normal_x"] = np.concatenate(invar["normal_x"], axis=0)
                 invar["normal_y"] = np.concatenate(invar["normal_y"], axis=0)
                 invar["normal_z"] = np.concatenate(invar["normal_z"], axis=0)
-                invar["area"] = np.concatenate(invar["area"], axis=0)
+                # Compute area from the original mesh
+                invar["area"] = np.ones_like(invar["x"]) * (
+                    np.sum(triangle_areas) / nr_points
+                )
 
                 # sample from the param ranges
                 params = parameterization.sample(nr_points, quasirandom=quasirandom)
@@ -234,29 +235,14 @@ def _sample_triangle(
 
 
 # area of array of triangles
-def _area_of_triangles(
-    v0, v1, v2
-):  # ref https://math.stackexchange.com/questions/128991/how-to-calculate-the-area-of-a-3d-triangle
-    a = np.sqrt(
-        (v0[:, 0] - v1[:, 0]) ** 2
-        + (v0[:, 1] - v1[:, 1]) ** 2
-        + (v0[:, 2] - v1[:, 2]) ** 2
-        + 1e-10
-    )
-    b = np.sqrt(
-        (v1[:, 0] - v2[:, 0]) ** 2
-        + (v1[:, 1] - v2[:, 1]) ** 2
-        + (v1[:, 2] - v2[:, 2]) ** 2
-        + 1e-10
-    )
-    c = np.sqrt(
-        (v0[:, 0] - v2[:, 0]) ** 2
-        + (v0[:, 1] - v2[:, 1]) ** 2
-        + (v0[:, 2] - v2[:, 2]) ** 2
-        + 1e-10
-    )
-    s = (a + b + c) / 2
-    area = np.sqrt(s * (s - a) * (s - b) * (s - c) + 1e-10)
+def _area_of_triangles(v0, v1, v2):
+    edge1 = v1 - v0
+    edge2 = v2 - v0
+
+    # use cross product to find the area of the triangle
+    cross_product = np.cross(edge1, edge2)
+    area = np.linalg.norm(cross_product, axis=1) / 2
+
     return area
 
 
diff --git a/test/test_pdes/test_linear_elasticity.py b/test/test_pdes/test_linear_elasticity.py
index fac5fe47..c92e20c5 100644
--- a/test/test_pdes/test_linear_elasticity.py
+++ b/test/test_pdes/test_linear_elasticity.py
@@ -294,9 +294,15 @@ def test_linear_elasticity_equations():
     assert np.allclose(traction_x_eval_pred, traction_x_true), "Test Failed!"
     assert np.allclose(traction_y_eval_pred, traction_y_true), "Test Failed!"
     assert np.allclose(traction_z_eval_pred, traction_z_true), "Test Failed!"
-    assert np.allclose(navier_x_eval_pred, navier_x_true, rtol=1e-3), "Test Failed!"
-    assert np.allclose(navier_y_eval_pred, navier_y_true, rtol=1e-3), "Test Failed!"
-    assert np.allclose(navier_z_eval_pred, navier_z_true, rtol=1e-3), "Test Failed!"
+    assert np.allclose(
+        navier_x_eval_pred, navier_x_true, atol=1e-3, rtol=1e-2
+    ), "Test Failed!"
+    assert np.allclose(
+        navier_y_eval_pred, navier_y_true, atol=1e-3, rtol=1e-2
+    ), "Test Failed!"
+    assert np.allclose(
+        navier_z_eval_pred, navier_z_true, atol=1e-3, rtol=1e-2
+    ), "Test Failed!"
 
 
 def test_linear_elasticity_plane_stress_equations():