From 952a6dcd5dfaa8d8cec794b6f3e322226cd69b8d Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 16 Jun 2021 23:25:03 -0500 Subject: [PATCH] Reproduce two minima in _make_cross_face_batch Gauss-Newton --- .../connection/opposite_face.py | 51 +++++++++++++++++-- 1 file changed, 47 insertions(+), 4 deletions(-) diff --git a/meshmode/discretization/connection/opposite_face.py b/meshmode/discretization/connection/opposite_face.py index b3f2d3a07..5b7849734 100644 --- a/meshmode/discretization/connection/opposite_face.py +++ b/meshmode/discretization/connection/opposite_face.py @@ -67,10 +67,13 @@ def _make_cross_face_batches(actx, src_mesh_grp = src_bdry_discr.mesh.groups[i_src_grp] src_grp = src_bdry_discr.groups[i_src_grp] - src_unit_nodes = _find_src_unit_nodes_by_matching( - tgt_bdry_nodes=tgt_bdry_nodes, - src_bdry_nodes=src_bdry_nodes, - src_grp=src_grp, tol=tol) + if 0: + src_unit_nodes = _find_src_unit_nodes_by_matching( + tgt_bdry_nodes=tgt_bdry_nodes, + src_bdry_nodes=src_bdry_nodes, + src_grp=src_grp, tol=tol) + else: + src_unit_nodes = None if src_unit_nodes is None: src_unit_nodes = _find_src_unit_nodes_via_gauss_newton( tgt_bdry_nodes=tgt_bdry_nodes, @@ -283,6 +286,46 @@ def get_map_jacobian(unit_nodes): max_resid = np.max(np.abs(resid)) if max_resid < tol: + if 1: + near_zero = np.abs(src_unit_nodes) < 1e-5 + near_one = np.abs(src_unit_nodes-1) < 1e-5 + near_mone = np.abs(src_unit_nodes+1) < 1e-5 + reasonable = near_zero | near_one | near_mone + if np.any(~reasonable): + unreasonable_els, unreasonable_nodes = np.where( + np.any(~reasonable, axis=0)) + mapped = apply_map(src_unit_nodes) + # print(src_unit_nodes[~reasonable]) + urel = unreasonable_els[0] + tgt_urel = tgt_bdry_nodes[:, urel] + src_urel = src_bdry_nodes[:, urel] + + dist_vecs = ( + tgt_urel.reshape(3, -1, 1) - src_urel.reshape(3, 1, -1)) + dists = np.sqrt(np.sum(dist_vecs**2, axis=0)) + + mat = (dists < 1e-12) + from_indices = np.array([np.where(row)[0][0] for row in mat]) + + sunodes = src_bdry_discr.groups[src_grp.index].unit_nodes + + print("GOOD UNIT NODES") + print(sunodes[:, from_indices]) + print("BAD UNIT NODES") + print(src_unit_nodes[:, urel]) + print() + + fixed_src_unit_nodes = src_unit_nodes.copy() + fixed_src_unit_nodes[:, urel] = sunodes[:, from_indices] + mapped_fixed = apply_map(fixed_src_unit_nodes) + + print("FIXED RESIDUAL") + print(np.max(np.abs(mapped_fixed[:, urel] + - tgt_bdry_nodes[:, urel]))) + print("UNFIXED RESIDUAL") + print(np.max(np.abs(mapped[:, urel] - tgt_bdry_nodes[:, urel]))) + 1/0 + logger.debug("_find_src_unit_nodes_via_gauss_newton: done, " "final residual: %g", max_resid) return src_unit_nodes