From cb97c26a29a6acb8559b6dbeb7316e8a2b145b2e Mon Sep 17 00:00:00 2001
From: Matt Wala <wala1@illinois.edu>
Date: Tue, 20 Sep 2016 14:32:11 -0500
Subject: [PATCH] check_nodal_adj_against_geometry(): Use least squares to find
 the barycentric coordinates.

---
 meshmode/mesh/refinement/utils.py | 13 ++++++++++---
 test/test_refinement.py           |  4 ++--
 2 files changed, 12 insertions(+), 5 deletions(-)

diff --git a/meshmode/mesh/refinement/utils.py b/meshmode/mesh/refinement/utils.py
index d34fa14a..936f9107 100644
--- a/meshmode/mesh/refinement/utils.py
+++ b/meshmode/mesh/refinement/utils.py
@@ -80,7 +80,7 @@ def check_nodal_adj_against_geometry(mesh, tol=1e-12):
                     nearby_origin_vertex = mesh.vertices[
                             :, nearby_grp.vertex_indices[nearby_iel][0]]  # noqa
                     transformation = np.empty(
-                            (mesh.ambient_dim, mesh.ambient_dim))
+                            (len(mesh.vertices), nvertices_per_element-1))
                     vertex_transformed = vertex - nearby_origin_vertex
 
                     for inearby_vertex_index, nearby_vertex_index in enumerate(
@@ -88,10 +88,17 @@ def check_nodal_adj_against_geometry(mesh, tol=1e-12):
                         nearby_vertex = mesh.vertices[:, nearby_vertex_index]
                         transformation[:, inearby_vertex_index] = \
                                 nearby_vertex - nearby_origin_vertex
-                    bary_coord = np.linalg.solve(transformation, vertex_transformed)
+                    bary_coord, residual = \
+                            np.linalg.lstsq(transformation, vertex_transformed)[0:2]
+
+                    is_in_element_span = (
+                            residual.size == 0 or
+                            np.linalg.norm(vertex_transformed) == 0 or
+                            residual / np.linalg.norm(vertex_transformed) <= tol)
 
                     is_connected = (
-                            np.sum(bary_coord) <= 1+tol
+                            is_in_element_span
+                            and np.sum(bary_coord) <= 1+tol
                             and (bary_coord >= -tol).all())
                     el1 = group_and_iel_to_global_iel(nearby_igrp, nearby_iel)
                     el2 = group_and_iel_to_global_iel(igrp, iel_grp)
diff --git a/test/test_refinement.py b/test/test_refinement.py
index fca3aa93..bed61eb9 100644
--- a/test/test_refinement.py
+++ b/test/test_refinement.py
@@ -84,8 +84,8 @@ def uniform_refine_flags(mesh):
         partial(
             make_curve_mesh,
             partial(ellipse, 3),
-            np.linspace(0, 1, 20),
-            order=2),
+            np.linspace(0, 1, 21),
+            order=1),
         uniform_refine_flags,
         4),
 
-- 
GitLab