From 89a710ecba9f1c21e4fb08e7a3b3c7e934a3dc0e Mon Sep 17 00:00:00 2001 From: Jack Betteridge Date: Wed, 25 Sep 2024 17:23:21 +0100 Subject: [PATCH 01/20] WIP --- ngsPETSc/utils/firedrake/meshes.py | 83 ++++++++++++++++-------------- 1 file changed, 44 insertions(+), 39 deletions(-) diff --git a/ngsPETSc/utils/firedrake/meshes.py b/ngsPETSc/utils/firedrake/meshes.py index ff4716e..37c1515 100644 --- a/ngsPETSc/utils/firedrake/meshes.py +++ b/ngsPETSc/utils/firedrake/meshes.py @@ -78,51 +78,57 @@ def curveField(self, order, tol=1e-8): ''' #Checking if the mesh is a surface mesh or two dimensional mesh - surf = len(self.netgen_mesh.Elements3D()) == 0 + surf = len(self.netgen_mesh.Elements3D()) == 0 # REMOVE + if len(self.netgen_mesh.Elements3D()) == 0: + ng_element = self.netgen_mesh.Elements2D + else: + ng_element = self.netgen_mesh.Elements3D + ng_dimension = len(ng_element()) + #Constructing mesh as a function low_order_element = self.coordinates.function_space().ufl_element().sub_elements[0] - element = low_order_element.reconstruct(degree=order) - space = fd.VectorFunctionSpace(self, fd.BrokenElement(element)) - newFunctionCoordinates = fd.assemble(interpolate(self.coordinates, space)) + ufl_element = low_order_element.reconstruct(degree=order) + firedrake_space = fd.VectorFunctionSpace(self, fd.BrokenElement(ufl_element)) + newFunctionCoordinates = fd.assemble(interpolate(self.coordinates, firedrake_space)) + #Computing reference points using fiat fiat_element = newFunctionCoordinates.function_space().finat_element.fiat_equivalent entity_ids = fiat_element.entity_dofs() nodes = fiat_element.dual_basis() - refPts = [] + ref = [] for dim in entity_ids: for entity in entity_ids[dim]: for dof in entity_ids[dim][entity]: # Assert singleton point for each node. pt, = nodes[dof].get_point_dict().keys() - refPts.append(pt) - V = newFunctionCoordinates.dat.data - refPts = np.array(refPts) - els = {True: self.netgen_mesh.Elements2D, False: self.netgen_mesh.Elements3D} + ref.append(pt) + reference_space_points = np.array(ref) + #Mapping to the physical domain + els = {True: self.netgen_mesh.Elements2D, False: self.netgen_mesh.Elements3D} # REMOVE + physical_space_points = np.ndarray((ng_dimension, reference_space_points.shape[0], self.geometric_dimension())) + curved_space_points = np.ndarray((ng_dimension, reference_space_points.shape[0], self.geometric_dimension())) + if self.comm.rank == 0: - physPts = np.ndarray((len(els[surf]()), - refPts.shape[0], self.geometric_dimension())) - self.netgen_mesh.CalcElementMapping(refPts, physPts) - #Cruving the mesh + #Curving the mesh on rank 0 + self.netgen_mesh.CalcElementMapping(reference_space_points, physical_space_points) self.netgen_mesh.Curve(order) - curvedPhysPts = np.ndarray((len(els[surf]()), - refPts.shape[0], self.geometric_dimension())) - self.netgen_mesh.CalcElementMapping(refPts, curvedPhysPts) - curved = els[surf]().NumPy()["curved"] + self.netgen_mesh.CalcElementMapping(reference_space_points, curved_space_points) + curved = ng_element().NumPy()["curved"] else: - physPts = np.ndarray((len(els[surf]()), - refPts.shape[0], self.geometric_dimension())) - curvedPhysPts = np.ndarray((len(els[surf]()), - refPts.shape[0], self.geometric_dimension())) - curved = np.array((len(els[surf]()),1)) - physPts = self.comm.bcast(physPts, root=0) - curvedPhysPts = self.comm.bcast(curvedPhysPts, root=0) + curved = np.array((ng_dimension, 1)) + + # Broadcast curving + physical_space_points = self.comm.bcast(physical_space_points, root=0) + curved_space_points = self.comm.bcast(curved_space_points, root=0) curved = self.comm.bcast(curved, root=0) cellMap = newFunctionCoordinates.cell_node_map() - for i in range(physPts.shape[0]): + + + for i in range(physical_space_points.shape[0]): #Inefficent code but runs only on curved elements if curved[i]: - pts = physPts[i][0:refPts.shape[0]] + pts = physical_space_points[i][0:reference_space_points.shape[0]] bary = sum([np.array(pts[i]) for i in range(len(pts))])/len(pts) Idx = self.locate_cell(bary) isInMesh = (0<=Idx tol: - fd.logging.warning("[{}] Not able to curve Firedrake element {} \ - ({}) -- residual: {}".format(self.comm.rank, Idx,i, res)) + fd.logging.warning(f"[{self.comm.rank}] Not able to curve Firedrake element {Idx} ({i}) -- residual: {res}") else: - for j, datIdx in enumerate(cellMap.values[Idx][0:refPts.shape[0]]): + for j, datIdx in enumerate(cellMap.values[Idx][0:reference_space_points.shape[0]]): for dim in range(self.geometric_dimension()): - coo = curvedPhysPts[i][j][dim] + coo = curved_space_points[i][j][dim] newFunctionCoordinates.sub(dim).dat.data[datIdx] = coo else: if isInMesh: p = [np.argmin(np.sum((pts - pt)**2, axis=1)) - for pt in V[cellMap.values[Idx]][0:refPts.shape[0]]] - curvedPhysPts[i] = curvedPhysPts[i][p] - res = np.linalg.norm(pts[p]-V[cellMap.values[Idx]][0:refPts.shape[0]]) + for pt in newFunctionCoordinates.dat.data[cellMap.values[Idx]][0:reference_space_points.shape[0]]] + curved_space_points[i] = curved_space_points[i][p] + res = np.linalg.norm(pts[p]-newFunctionCoordinates.dat.data[cellMap.values[Idx]][0:reference_space_points.shape[0]]) else: res = np.inf res = self.comm.gather(res, root=0) @@ -160,11 +165,11 @@ def curveField(self, order, tol=1e-8): fd.logging.warning("[{}, {}] Not able to curve Firedrake element {} \ ({}) -- residual: {}".format(self.comm.rank, shared, Idx,i, res)) else: - for j, datIdx in enumerate(cellMap.values[Idx][0:refPts.shape[0]]): + for j, datIdx in enumerate(cellMap.values[Idx][0:reference_space_points.shape[0]]): for dim in range(self.geometric_dimension()): - coo = curvedPhysPts[i][j][dim] + coo = curved_space_points[i][j][dim] newFunctionCoordinates.sub(dim).dat.data[datIdx] = coo - + # ~ breakpoint() return newFunctionCoordinates def splitToQuads(plex, dim, comm): From 814ff669dbc9feb8f64d37af08b2716230d19936 Mon Sep 17 00:00:00 2001 From: Jack Betteridge Date: Thu, 26 Sep 2024 01:05:55 +0100 Subject: [PATCH 02/20] WIP: almost loop free --- ngsPETSc/utils/firedrake/meshes.py | 120 ++++++++++++++++++----------- 1 file changed, 74 insertions(+), 46 deletions(-) diff --git a/ngsPETSc/utils/firedrake/meshes.py b/ngsPETSc/utils/firedrake/meshes.py index 37c1515..9a5a605 100644 --- a/ngsPETSc/utils/firedrake/meshes.py +++ b/ngsPETSc/utils/firedrake/meshes.py @@ -118,57 +118,85 @@ def curveField(self, order, tol=1e-8): else: curved = np.array((ng_dimension, 1)) - # Broadcast curving + # Broadcast curving data physical_space_points = self.comm.bcast(physical_space_points, root=0) curved_space_points = self.comm.bcast(curved_space_points, root=0) curved = self.comm.bcast(curved, root=0) - cellMap = newFunctionCoordinates.cell_node_map() + cell_node_map = newFunctionCoordinates.cell_node_map() + # Select only the points in curved cells + physical_space_points = physical_space_points[curved] + curved_space_points = curved_space_points[curved] + barycentres = np.average(physical_space_points, axis=1) + ng_index = [*map(self.locate_cell, barycentres)] + owned = [(0 <= ii < len(cell_node_map.values)) if ii is not None else False for ii in ng_index] - for i in range(physical_space_points.shape[0]): - #Inefficent code but runs only on curved elements - if curved[i]: - pts = physical_space_points[i][0:reference_space_points.shape[0]] - bary = sum([np.array(pts[i]) for i in range(len(pts))])/len(pts) - Idx = self.locate_cell(bary) - isInMesh = (0<=Idx tol: - fd.logging.warning(f"[{self.comm.rank}] Not able to curve Firedrake element {Idx} ({i}) -- residual: {res}") - else: - for j, datIdx in enumerate(cellMap.values[Idx][0:reference_space_points.shape[0]]): - for dim in range(self.geometric_dimension()): - coo = curved_space_points[i][j][dim] - newFunctionCoordinates.sub(dim).dat.data[datIdx] = coo - else: - if isInMesh: - p = [np.argmin(np.sum((pts - pt)**2, axis=1)) - for pt in newFunctionCoordinates.dat.data[cellMap.values[Idx]][0:reference_space_points.shape[0]]] - curved_space_points[i] = curved_space_points[i][p] - res = np.linalg.norm(pts[p]-newFunctionCoordinates.dat.data[cellMap.values[Idx]][0:reference_space_points.shape[0]]) - else: - res = np.inf - res = self.comm.gather(res, root=0) - res = self.comm.bcast(res, root=0) - rank = np.argmin(res) - if self.comm.rank == rank: - if res[rank] > tol: - fd.logging.warning("[{}, {}] Not able to curve Firedrake element {} \ - ({}) -- residual: {}".format(self.comm.rank, shared, Idx,i, res)) - else: - for j, datIdx in enumerate(cellMap.values[Idx][0:reference_space_points.shape[0]]): - for dim in range(self.geometric_dimension()): - coo = curved_space_points[i][j][dim] - newFunctionCoordinates.sub(dim).dat.data[datIdx] = coo + # Select only the points owned by this rank + physical_space_points = physical_space_points[owned] + curved_space_points = curved_space_points[owned] + barycentres = barycentres[owned] + ng_index = [idx for idx, o in zip(ng_index, owned) if o] + + # Do we want the min or max of these??? + norms = np.linalg.norm(physical_space_points - newFunctionCoordinates.dat.data[cell_node_map.values[ng_index]], axis=2) + + # PyOP2 index + pyop2_index = [] + for ngidx in ng_index: + pyop2_index.extend(cell_node_map.values[ngidx]) + np.array(pyop2_index) + + for dim in range(self.geometric_dimension()): + newFunctionCoordinates.sub(dim).dat.data[pyop2_index] = curved_space_points[:,:,dim].flatten() + + # ~ breakpoint() + + # ~ for i in range(physical_space_points.shape[0]): + # ~ #Inefficent code but runs only on curved elements + # ~ if curved[i]: + # ~ pts = physical_space_points[i][0:reference_space_points.shape[0]] + # ~ bary = sum([np.array(pts[i]) for i in range(len(pts))])/len(pts) + # ~ Idx = self.locate_cell(bary) + # ~ isInMesh = (0<=Idx tol: + # ~ fd.logging.warning(f"[{self.comm.rank}] Not able to curve Firedrake element {Idx} ({i}) -- residual: {res}") + # ~ else: + # ~ for j, datIdx in enumerate(cell_node_map.values[Idx][0:reference_space_points.shape[0]]): + # ~ for dim in range(self.geometric_dimension()): + # ~ coo = curved_space_points[i][j][dim] + # ~ newFunctionCoordinates.sub(dim).dat.data[datIdx] = coo + # ~ else: + # ~ if isInMesh: + # ~ p = [np.argmin(np.sum((pts - pt)**2, axis=1)) + # ~ for pt in newFunctionCoordinates.dat.data[cell_node_map.values[Idx]][0:reference_space_points.shape[0]]] + # ~ curved_space_points[i] = curved_space_points[i][p] + # ~ res = np.linalg.norm(pts[p]-newFunctionCoordinates.dat.data[cell_node_map.values[Idx]][0:reference_space_points.shape[0]]) + # ~ else: + # ~ res = np.inf + # ~ res = self.comm.gather(res, root=0) + # ~ res = self.comm.bcast(res, root=0) + # ~ rank = np.argmin(res) + # ~ if self.comm.rank == rank: + # ~ if res[rank] > tol: + # ~ fd.logging.warning("[{}, {}] Not able to curve Firedrake element {} \ + # ~ ({}) -- residual: {}".format(self.comm.rank, shared, Idx,i, res)) + # ~ else: + # ~ for j, datIdx in enumerate(cell_node_map.values[Idx][0:reference_space_points.shape[0]]): + # ~ for dim in range(self.geometric_dimension()): + # ~ coo = curved_space_points[i][j][dim] + # ~ newFunctionCoordinates.sub(dim).dat.data[datIdx] = coo # ~ breakpoint() return newFunctionCoordinates From 4bb63af0e97b353ceee6d1eb091c27da5f849e4a Mon Sep 17 00:00:00 2001 From: Jack Betteridge Date: Tue, 1 Oct 2024 17:45:37 +0100 Subject: [PATCH 03/20] Additional tidying, difficult numpy indexing rules --- ngsPETSc/utils/firedrake/meshes.py | 128 +++++++++++++++-------------- 1 file changed, 66 insertions(+), 62 deletions(-) diff --git a/ngsPETSc/utils/firedrake/meshes.py b/ngsPETSc/utils/firedrake/meshes.py index 9a5a605..5a6e5ff 100644 --- a/ngsPETSc/utils/firedrake/meshes.py +++ b/ngsPETSc/utils/firedrake/meshes.py @@ -12,6 +12,12 @@ import numpy as np from petsc4py import PETSc +try: + from scipy.spatial.distance import cdist + HAVE_SCIPY = True +except ImportError: + HAVE_SCIPY = False + import netgen import netgen.meshing as ngm from netgen.meshing import MeshingParameters @@ -70,6 +76,45 @@ def refineMarkedElements(self, mark): else: raise NotImplementedError("No implementation for dimension other than 2 and 3.") + +def _slow_cdist(XA, XB): + dist = np.zeros([len(XA), len(XB)]) + for ii, a in enumerate(XA): + for jj, b in enumerate(XB): + dist[ii, jj] = np.linalg.norm(b - a) + return dist + + +if not HAVE_SCIPY: + cdist = _slow_cdist + + +def find_permutation(points_a, points_b, tol=1e-5): + """ Find all permutations between a list of two sets of points. + + Given two numpy arrays of shape (ncells, npoints, dim) containing + floating point coordinates for each cell, determine each index + permutation that takes `points_a` to `points_b`. Ie: + ``` + permutation = find_permutation(points_a, points_b) + assert np.allclose(points_a[permutation], points_b, rtol=0, atol=tol) + ``` + """ + if points_a.shape != points_b.shape: + raise ValueError("`points_a` and `points_b` must have the same shape.") + + p = [np.where(cdist(a, b).T < tol)[1] for a, b in zip(points_a, points_b)] + try: + permutation = np.array(p, ndmin=2) + except ValueError: + raise ValueError("It was not possible to find a permutation for every cell within the provided tolerance") + + if permutation.shape != points_a.shape[0:2]: + raise ValueError("It was not possible to find a permutation for every cell within the provided tolerance") + + return permutation + + def curveField(self, order, tol=1e-8): ''' This method returns a curved mesh as a Firedrake function. @@ -78,21 +123,21 @@ def curveField(self, order, tol=1e-8): ''' #Checking if the mesh is a surface mesh or two dimensional mesh - surf = len(self.netgen_mesh.Elements3D()) == 0 # REMOVE if len(self.netgen_mesh.Elements3D()) == 0: ng_element = self.netgen_mesh.Elements2D else: ng_element = self.netgen_mesh.Elements3D ng_dimension = len(ng_element()) + geom_dim = self.geometric_dimension() #Constructing mesh as a function low_order_element = self.coordinates.function_space().ufl_element().sub_elements[0] ufl_element = low_order_element.reconstruct(degree=order) firedrake_space = fd.VectorFunctionSpace(self, fd.BrokenElement(ufl_element)) - newFunctionCoordinates = fd.assemble(interpolate(self.coordinates, firedrake_space)) + new_coordinates = fd.assemble(interpolate(self.coordinates, firedrake_space)) #Computing reference points using fiat - fiat_element = newFunctionCoordinates.function_space().finat_element.fiat_equivalent + fiat_element = new_coordinates.function_space().finat_element.fiat_equivalent entity_ids = fiat_element.entity_dofs() nodes = fiat_element.dual_basis() ref = [] @@ -105,9 +150,8 @@ def curveField(self, order, tol=1e-8): reference_space_points = np.array(ref) #Mapping to the physical domain - els = {True: self.netgen_mesh.Elements2D, False: self.netgen_mesh.Elements3D} # REMOVE - physical_space_points = np.ndarray((ng_dimension, reference_space_points.shape[0], self.geometric_dimension())) - curved_space_points = np.ndarray((ng_dimension, reference_space_points.shape[0], self.geometric_dimension())) + physical_space_points = np.ndarray((ng_dimension, reference_space_points.shape[0], geom_dim)) + curved_space_points = np.ndarray((ng_dimension, reference_space_points.shape[0], geom_dim)) if self.comm.rank == 0: #Curving the mesh on rank 0 @@ -122,7 +166,7 @@ def curveField(self, order, tol=1e-8): physical_space_points = self.comm.bcast(physical_space_points, root=0) curved_space_points = self.comm.bcast(curved_space_points, root=0) curved = self.comm.bcast(curved, root=0) - cell_node_map = newFunctionCoordinates.cell_node_map() + cell_node_map = new_coordinates.cell_node_map() # Select only the points in curved cells physical_space_points = physical_space_points[curved] @@ -137,8 +181,7 @@ def curveField(self, order, tol=1e-8): barycentres = barycentres[owned] ng_index = [idx for idx, o in zip(ng_index, owned) if o] - # Do we want the min or max of these??? - norms = np.linalg.norm(physical_space_points - newFunctionCoordinates.dat.data[cell_node_map.values[ng_index]], axis=2) + breakpoint() # PyOP2 index pyop2_index = [] @@ -146,59 +189,20 @@ def curveField(self, order, tol=1e-8): pyop2_index.extend(cell_node_map.values[ngidx]) np.array(pyop2_index) - for dim in range(self.geometric_dimension()): - newFunctionCoordinates.sub(dim).dat.data[pyop2_index] = curved_space_points[:,:,dim].flatten() - - # ~ breakpoint() - - # ~ for i in range(physical_space_points.shape[0]): - # ~ #Inefficent code but runs only on curved elements - # ~ if curved[i]: - # ~ pts = physical_space_points[i][0:reference_space_points.shape[0]] - # ~ bary = sum([np.array(pts[i]) for i in range(len(pts))])/len(pts) - # ~ Idx = self.locate_cell(bary) - # ~ isInMesh = (0<=Idx tol: - # ~ fd.logging.warning(f"[{self.comm.rank}] Not able to curve Firedrake element {Idx} ({i}) -- residual: {res}") - # ~ else: - # ~ for j, datIdx in enumerate(cell_node_map.values[Idx][0:reference_space_points.shape[0]]): - # ~ for dim in range(self.geometric_dimension()): - # ~ coo = curved_space_points[i][j][dim] - # ~ newFunctionCoordinates.sub(dim).dat.data[datIdx] = coo - # ~ else: - # ~ if isInMesh: - # ~ p = [np.argmin(np.sum((pts - pt)**2, axis=1)) - # ~ for pt in newFunctionCoordinates.dat.data[cell_node_map.values[Idx]][0:reference_space_points.shape[0]]] - # ~ curved_space_points[i] = curved_space_points[i][p] - # ~ res = np.linalg.norm(pts[p]-newFunctionCoordinates.dat.data[cell_node_map.values[Idx]][0:reference_space_points.shape[0]]) - # ~ else: - # ~ res = np.inf - # ~ res = self.comm.gather(res, root=0) - # ~ res = self.comm.bcast(res, root=0) - # ~ rank = np.argmin(res) - # ~ if self.comm.rank == rank: - # ~ if res[rank] > tol: - # ~ fd.logging.warning("[{}, {}] Not able to curve Firedrake element {} \ - # ~ ({}) -- residual: {}".format(self.comm.rank, shared, Idx,i, res)) - # ~ else: - # ~ for j, datIdx in enumerate(cell_node_map.values[Idx][0:reference_space_points.shape[0]]): - # ~ for dim in range(self.geometric_dimension()): - # ~ coo = curved_space_points[i][j][dim] - # ~ newFunctionCoordinates.sub(dim).dat.data[datIdx] = coo - # ~ breakpoint() - return newFunctionCoordinates + # Find the correct coordinate permutation for each cell + permutation = find_permutation( + physical_space_points, + new_coordinates.dat.data[pyop2_index].reshape(physical_space_points.shape) + ) + + # Apply the permutation to each cell in turn + for ii, p in enumerate(curved_space_points): + curved_space_points[ii] = p[permutation[ii]] + + # Assign the curved coordinates to the dat + new_coordinates.dat.data[pyop2_index] = curved_space_points.reshape(-1, geom_dim) + + return new_coordinates def splitToQuads(plex, dim, comm): ''' From 0f27605135c9f8935f6d6d29f9c47047fb5b8576 Mon Sep 17 00:00:00 2001 From: Jack Betteridge Date: Tue, 1 Oct 2024 17:46:02 +0100 Subject: [PATCH 04/20] Left breakpoint behind --- ngsPETSc/utils/firedrake/meshes.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/ngsPETSc/utils/firedrake/meshes.py b/ngsPETSc/utils/firedrake/meshes.py index 5a6e5ff..50f11d9 100644 --- a/ngsPETSc/utils/firedrake/meshes.py +++ b/ngsPETSc/utils/firedrake/meshes.py @@ -181,8 +181,6 @@ def curveField(self, order, tol=1e-8): barycentres = barycentres[owned] ng_index = [idx for idx, o in zip(ng_index, owned) if o] - breakpoint() - # PyOP2 index pyop2_index = [] for ngidx in ng_index: From 3b1e1b53b5596b3ac1d57fae0a33d881da1e90ac Mon Sep 17 00:00:00 2001 From: Jack Betteridge Date: Tue, 1 Oct 2024 17:50:06 +0100 Subject: [PATCH 05/20] Lint --- ngsPETSc/utils/firedrake/meshes.py | 25 ++++++++++++++++--------- 1 file changed, 16 insertions(+), 9 deletions(-) diff --git a/ngsPETSc/utils/firedrake/meshes.py b/ngsPETSc/utils/firedrake/meshes.py index 50f11d9..32b6858 100644 --- a/ngsPETSc/utils/firedrake/meshes.py +++ b/ngsPETSc/utils/firedrake/meshes.py @@ -78,11 +78,11 @@ def refineMarkedElements(self, mark): def _slow_cdist(XA, XB): - dist = np.zeros([len(XA), len(XB)]) - for ii, a in enumerate(XA): - for jj, b in enumerate(XB): - dist[ii, jj] = np.linalg.norm(b - a) - return dist + dist = np.zeros([len(XA), len(XB)]) + for ii, a in enumerate(XA): + for jj, b in enumerate(XB): + dist[ii, jj] = np.linalg.norm(b - a) + return dist if not HAVE_SCIPY: @@ -106,11 +106,17 @@ def find_permutation(points_a, points_b, tol=1e-5): p = [np.where(cdist(a, b).T < tol)[1] for a, b in zip(points_a, points_b)] try: permutation = np.array(p, ndmin=2) - except ValueError: - raise ValueError("It was not possible to find a permutation for every cell within the provided tolerance") + except ValueError as e: + raise ValueError( + "It was not possible to find a permutation for every cell" + " within the provided tolerance" + ) from e if permutation.shape != points_a.shape[0:2]: - raise ValueError("It was not possible to find a permutation for every cell within the provided tolerance") + raise ValueError( + "It was not possible to find a permutation for every cell" + " within the provided tolerance" + ) return permutation @@ -190,7 +196,8 @@ def curveField(self, order, tol=1e-8): # Find the correct coordinate permutation for each cell permutation = find_permutation( physical_space_points, - new_coordinates.dat.data[pyop2_index].reshape(physical_space_points.shape) + new_coordinates.dat.data[pyop2_index].reshape(physical_space_points.shape), + tol=tol ) # Apply the permutation to each cell in turn From f3ad982ad83654284465ea53dc67928c90d0ae52 Mon Sep 17 00:00:00 2001 From: Jack Betteridge Date: Tue, 1 Oct 2024 18:18:04 +0100 Subject: [PATCH 06/20] More notes --- ngsPETSc/utils/firedrake/meshes.py | 23 +++++++++++++++-------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/ngsPETSc/utils/firedrake/meshes.py b/ngsPETSc/utils/firedrake/meshes.py index 32b6858..e8ee7fa 100644 --- a/ngsPETSc/utils/firedrake/meshes.py +++ b/ngsPETSc/utils/firedrake/meshes.py @@ -128,7 +128,7 @@ def curveField(self, order, tol=1e-8): :arg order: the order of the curved mesh ''' - #Checking if the mesh is a surface mesh or two dimensional mesh + # Check if the mesh is a surface mesh or two dimensional mesh if len(self.netgen_mesh.Elements3D()) == 0: ng_element = self.netgen_mesh.Elements2D else: @@ -136,13 +136,13 @@ def curveField(self, order, tol=1e-8): ng_dimension = len(ng_element()) geom_dim = self.geometric_dimension() - #Constructing mesh as a function + # Construct the mesh as a Firedrake function low_order_element = self.coordinates.function_space().ufl_element().sub_elements[0] ufl_element = low_order_element.reconstruct(degree=order) firedrake_space = fd.VectorFunctionSpace(self, fd.BrokenElement(ufl_element)) new_coordinates = fd.assemble(interpolate(self.coordinates, firedrake_space)) - #Computing reference points using fiat + # Compute reference points using fiat fiat_element = new_coordinates.function_space().finat_element.fiat_equivalent entity_ids = fiat_element.entity_dofs() nodes = fiat_element.dual_basis() @@ -155,12 +155,17 @@ def curveField(self, order, tol=1e-8): ref.append(pt) reference_space_points = np.array(ref) - #Mapping to the physical domain - physical_space_points = np.ndarray((ng_dimension, reference_space_points.shape[0], geom_dim)) - curved_space_points = np.ndarray((ng_dimension, reference_space_points.shape[0], geom_dim)) + # Map to the physical domain + physical_space_points = np.ndarray( + (ng_dimension, reference_space_points.shape[0], geom_dim) + ) + curved_space_points = np.ndarray( + (ng_dimension, reference_space_points.shape[0], geom_dim) + ) + # Curve the mesh on rank 0 only + # JBTODO: (why?) if self.comm.rank == 0: - #Curving the mesh on rank 0 self.netgen_mesh.CalcElementMapping(reference_space_points, physical_space_points) self.netgen_mesh.Curve(order) self.netgen_mesh.CalcElementMapping(reference_space_points, curved_space_points) @@ -187,13 +192,15 @@ def curveField(self, order, tol=1e-8): barycentres = barycentres[owned] ng_index = [idx for idx, o in zip(ng_index, owned) if o] - # PyOP2 index + # Get the PyOP2 indices corresponding to the netgen indices pyop2_index = [] for ngidx in ng_index: pyop2_index.extend(cell_node_map.values[ngidx]) np.array(pyop2_index) # Find the correct coordinate permutation for each cell + # JBTODO: This should be moved to the next loop if we have to loop + # over cells any way to actually perform the permutation permutation = find_permutation( physical_space_points, new_coordinates.dat.data[pyop2_index].reshape(physical_space_points.shape), From 37e96d913c114445d60b48c1d1a60968da9ff689 Mon Sep 17 00:00:00 2001 From: Jack Betteridge Date: Thu, 3 Oct 2024 01:47:17 +0100 Subject: [PATCH 07/20] This is a mess, revert and try again --- ngsPETSc/utils/firedrake/meshes.py | 106 ++++++++++++++++------------- 1 file changed, 57 insertions(+), 49 deletions(-) diff --git a/ngsPETSc/utils/firedrake/meshes.py b/ngsPETSc/utils/firedrake/meshes.py index e8ee7fa..3ff86e8 100644 --- a/ngsPETSc/utils/firedrake/meshes.py +++ b/ngsPETSc/utils/firedrake/meshes.py @@ -90,10 +90,10 @@ def _slow_cdist(XA, XB): def find_permutation(points_a, points_b, tol=1e-5): - """ Find all permutations between a list of two sets of points. + """ Find permutation between two sets of points. - Given two numpy arrays of shape (ncells, npoints, dim) containing - floating point coordinates for each cell, determine each index + Given two numpy arrays of shape (npoints, dim) containing + floating point coordinates for a cell, determine each index permutation that takes `points_a` to `points_b`. Ie: ``` permutation = find_permutation(points_a, points_b) @@ -103,18 +103,11 @@ def find_permutation(points_a, points_b, tol=1e-5): if points_a.shape != points_b.shape: raise ValueError("`points_a` and `points_b` must have the same shape.") - p = [np.where(cdist(a, b).T < tol)[1] for a, b in zip(points_a, points_b)] - try: - permutation = np.array(p, ndmin=2) - except ValueError as e: - raise ValueError( - "It was not possible to find a permutation for every cell" - " within the provided tolerance" - ) from e + permutation = np.where(cdist(points_a, points_b).T < tol)[1] - if permutation.shape != points_a.shape[0:2]: + if permutation.shape[0] != points_a.shape[0]: raise ValueError( - "It was not possible to find a permutation for every cell" + "It was not possible to find a permutation for the cell" " within the provided tolerance" ) @@ -155,64 +148,79 @@ def curveField(self, order, tol=1e-8): ref.append(pt) reference_space_points = np.array(ref) - # Map to the physical domain - physical_space_points = np.ndarray( - (ng_dimension, reference_space_points.shape[0], geom_dim) - ) - curved_space_points = np.ndarray( - (ng_dimension, reference_space_points.shape[0], geom_dim) - ) - # Curve the mesh on rank 0 only - # JBTODO: (why?) if self.comm.rank == 0: + # Construct numpy arrays for physical domain data + physical_space_points = np.ndarray( + (ng_dimension, reference_space_points.shape[0], geom_dim) + ) + curved_space_points = np.ndarray( + (ng_dimension, reference_space_points.shape[0], geom_dim) + ) self.netgen_mesh.CalcElementMapping(reference_space_points, physical_space_points) self.netgen_mesh.Curve(order) self.netgen_mesh.CalcElementMapping(reference_space_points, curved_space_points) curved = ng_element().NumPy()["curved"] + # Broadcast curving data + curved = self.comm.bcast(curved, root=0) + physical_space_points = physical_space_points[curved] + curved_space_points = curved_space_points[curved] else: - curved = np.array((ng_dimension, 1)) + curved = self.comm.bcast(None, root=0) + # Construct numpy arrays as buffers to receive physical domain data + ncurved = np.sum(curved) + physical_space_points = np.ndarray( + (ncurved, reference_space_points.shape[0], geom_dim) + ) + curved_space_points = np.ndarray( + (ncurved, reference_space_points.shape[0], geom_dim) + ) - # Broadcast curving data - physical_space_points = self.comm.bcast(physical_space_points, root=0) - curved_space_points = self.comm.bcast(curved_space_points, root=0) - curved = self.comm.bcast(curved, root=0) + # ~ import pytest; pytest.set_trace() + # ~ breakpoint() + # Broadcast physical domain data for curved points + self.comm.Bcast(physical_space_points, root=0) + self.comm.Bcast(curved_space_points, root=0) cell_node_map = new_coordinates.cell_node_map() # Select only the points in curved cells - physical_space_points = physical_space_points[curved] - curved_space_points = curved_space_points[curved] barycentres = np.average(physical_space_points, axis=1) ng_index = [*map(self.locate_cell, barycentres)] + + # Select only the indices of points owned by this rank owned = [(0 <= ii < len(cell_node_map.values)) if ii is not None else False for ii in ng_index] - # Select only the points owned by this rank - physical_space_points = physical_space_points[owned] - curved_space_points = curved_space_points[owned] - barycentres = barycentres[owned] - ng_index = [idx for idx, o in zip(ng_index, owned) if o] + owned_physical_space_points = physical_space_points[owned] + owned_curved_space_points = curved_space_points[owned] + owned_barycentres = barycentres[owned] + owned_index = [idx for idx, o in zip(ng_index, owned) if o] # Get the PyOP2 indices corresponding to the netgen indices pyop2_index = [] - for ngidx in ng_index: + for ngidx in owned_index: pyop2_index.extend(cell_node_map.values[ngidx]) - np.array(pyop2_index) - - # Find the correct coordinate permutation for each cell - # JBTODO: This should be moved to the next loop if we have to loop - # over cells any way to actually perform the permutation - permutation = find_permutation( - physical_space_points, - new_coordinates.dat.data[pyop2_index].reshape(physical_space_points.shape), - tol=tol - ) + pyop2_index = np.array(pyop2_index).reshape(-1, physical_space_points.shape[1]) + + # ~ breakpoint() + # Find the correct coordinate permutation for this cell and apply the + # permutation to each cell in turn. For owned cells this must not fail. + for ii, (physical, index) in enumerate(zip(owned_physical_space_points, pyop2_index)): + try: + permutation = find_permutation( + physical, + new_coordinates.dat.data[index].reshape(physical.shape), + tol=tol + ) + owned_curved_space_points[ii] = physical[permutation] + except ValueError: + fd.logging.warning( + f"[{self.comm.rank}] Not able to curve Firedrake element" + f" {owned_index[ii]}" + ) - # Apply the permutation to each cell in turn - for ii, p in enumerate(curved_space_points): - curved_space_points[ii] = p[permutation[ii]] # Assign the curved coordinates to the dat - new_coordinates.dat.data[pyop2_index] = curved_space_points.reshape(-1, geom_dim) + new_coordinates.dat.data[pyop2_index.flatten()] = owned_curved_space_points.reshape(-1, geom_dim) return new_coordinates From adbb85fa86b2f4e74b8f5996b9e1144342424ad8 Mon Sep 17 00:00:00 2001 From: Jack Betteridge Date: Thu, 3 Oct 2024 11:52:15 +0100 Subject: [PATCH 08/20] Revert "This is a mess, revert and try again" This reverts commit 37e96d913c114445d60b48c1d1a60968da9ff689. --- ngsPETSc/utils/firedrake/meshes.py | 106 +++++++++++++---------------- 1 file changed, 49 insertions(+), 57 deletions(-) diff --git a/ngsPETSc/utils/firedrake/meshes.py b/ngsPETSc/utils/firedrake/meshes.py index 3ff86e8..e8ee7fa 100644 --- a/ngsPETSc/utils/firedrake/meshes.py +++ b/ngsPETSc/utils/firedrake/meshes.py @@ -90,10 +90,10 @@ def _slow_cdist(XA, XB): def find_permutation(points_a, points_b, tol=1e-5): - """ Find permutation between two sets of points. + """ Find all permutations between a list of two sets of points. - Given two numpy arrays of shape (npoints, dim) containing - floating point coordinates for a cell, determine each index + Given two numpy arrays of shape (ncells, npoints, dim) containing + floating point coordinates for each cell, determine each index permutation that takes `points_a` to `points_b`. Ie: ``` permutation = find_permutation(points_a, points_b) @@ -103,11 +103,18 @@ def find_permutation(points_a, points_b, tol=1e-5): if points_a.shape != points_b.shape: raise ValueError("`points_a` and `points_b` must have the same shape.") - permutation = np.where(cdist(points_a, points_b).T < tol)[1] + p = [np.where(cdist(a, b).T < tol)[1] for a, b in zip(points_a, points_b)] + try: + permutation = np.array(p, ndmin=2) + except ValueError as e: + raise ValueError( + "It was not possible to find a permutation for every cell" + " within the provided tolerance" + ) from e - if permutation.shape[0] != points_a.shape[0]: + if permutation.shape != points_a.shape[0:2]: raise ValueError( - "It was not possible to find a permutation for the cell" + "It was not possible to find a permutation for every cell" " within the provided tolerance" ) @@ -148,79 +155,64 @@ def curveField(self, order, tol=1e-8): ref.append(pt) reference_space_points = np.array(ref) + # Map to the physical domain + physical_space_points = np.ndarray( + (ng_dimension, reference_space_points.shape[0], geom_dim) + ) + curved_space_points = np.ndarray( + (ng_dimension, reference_space_points.shape[0], geom_dim) + ) + # Curve the mesh on rank 0 only + # JBTODO: (why?) if self.comm.rank == 0: - # Construct numpy arrays for physical domain data - physical_space_points = np.ndarray( - (ng_dimension, reference_space_points.shape[0], geom_dim) - ) - curved_space_points = np.ndarray( - (ng_dimension, reference_space_points.shape[0], geom_dim) - ) self.netgen_mesh.CalcElementMapping(reference_space_points, physical_space_points) self.netgen_mesh.Curve(order) self.netgen_mesh.CalcElementMapping(reference_space_points, curved_space_points) curved = ng_element().NumPy()["curved"] - # Broadcast curving data - curved = self.comm.bcast(curved, root=0) - physical_space_points = physical_space_points[curved] - curved_space_points = curved_space_points[curved] else: - curved = self.comm.bcast(None, root=0) - # Construct numpy arrays as buffers to receive physical domain data - ncurved = np.sum(curved) - physical_space_points = np.ndarray( - (ncurved, reference_space_points.shape[0], geom_dim) - ) - curved_space_points = np.ndarray( - (ncurved, reference_space_points.shape[0], geom_dim) - ) + curved = np.array((ng_dimension, 1)) - # ~ import pytest; pytest.set_trace() - # ~ breakpoint() - # Broadcast physical domain data for curved points - self.comm.Bcast(physical_space_points, root=0) - self.comm.Bcast(curved_space_points, root=0) + # Broadcast curving data + physical_space_points = self.comm.bcast(physical_space_points, root=0) + curved_space_points = self.comm.bcast(curved_space_points, root=0) + curved = self.comm.bcast(curved, root=0) cell_node_map = new_coordinates.cell_node_map() # Select only the points in curved cells + physical_space_points = physical_space_points[curved] + curved_space_points = curved_space_points[curved] barycentres = np.average(physical_space_points, axis=1) ng_index = [*map(self.locate_cell, barycentres)] - - # Select only the indices of points owned by this rank owned = [(0 <= ii < len(cell_node_map.values)) if ii is not None else False for ii in ng_index] - owned_physical_space_points = physical_space_points[owned] - owned_curved_space_points = curved_space_points[owned] - owned_barycentres = barycentres[owned] - owned_index = [idx for idx, o in zip(ng_index, owned) if o] + # Select only the points owned by this rank + physical_space_points = physical_space_points[owned] + curved_space_points = curved_space_points[owned] + barycentres = barycentres[owned] + ng_index = [idx for idx, o in zip(ng_index, owned) if o] # Get the PyOP2 indices corresponding to the netgen indices pyop2_index = [] - for ngidx in owned_index: + for ngidx in ng_index: pyop2_index.extend(cell_node_map.values[ngidx]) - pyop2_index = np.array(pyop2_index).reshape(-1, physical_space_points.shape[1]) - - # ~ breakpoint() - # Find the correct coordinate permutation for this cell and apply the - # permutation to each cell in turn. For owned cells this must not fail. - for ii, (physical, index) in enumerate(zip(owned_physical_space_points, pyop2_index)): - try: - permutation = find_permutation( - physical, - new_coordinates.dat.data[index].reshape(physical.shape), - tol=tol - ) - owned_curved_space_points[ii] = physical[permutation] - except ValueError: - fd.logging.warning( - f"[{self.comm.rank}] Not able to curve Firedrake element" - f" {owned_index[ii]}" - ) + np.array(pyop2_index) + + # Find the correct coordinate permutation for each cell + # JBTODO: This should be moved to the next loop if we have to loop + # over cells any way to actually perform the permutation + permutation = find_permutation( + physical_space_points, + new_coordinates.dat.data[pyop2_index].reshape(physical_space_points.shape), + tol=tol + ) + # Apply the permutation to each cell in turn + for ii, p in enumerate(curved_space_points): + curved_space_points[ii] = p[permutation[ii]] # Assign the curved coordinates to the dat - new_coordinates.dat.data[pyop2_index.flatten()] = owned_curved_space_points.reshape(-1, geom_dim) + new_coordinates.dat.data[pyop2_index] = curved_space_points.reshape(-1, geom_dim) return new_coordinates From 9e1e9cd9200ff8e1ae4c916406304246bf918a01 Mon Sep 17 00:00:00 2001 From: Jack Betteridge Date: Thu, 3 Oct 2024 17:07:01 +0100 Subject: [PATCH 09/20] Tolerance needs to be dropped for cell location --- ngsPETSc/utils/firedrake/meshes.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ngsPETSc/utils/firedrake/meshes.py b/ngsPETSc/utils/firedrake/meshes.py index e8ee7fa..78801b2 100644 --- a/ngsPETSc/utils/firedrake/meshes.py +++ b/ngsPETSc/utils/firedrake/meshes.py @@ -183,7 +183,7 @@ def curveField(self, order, tol=1e-8): physical_space_points = physical_space_points[curved] curved_space_points = curved_space_points[curved] barycentres = np.average(physical_space_points, axis=1) - ng_index = [*map(self.locate_cell, barycentres)] + ng_index = [*map(lambda x: self.locate_cell(x, tolerance=0.01), barycentres)] owned = [(0 <= ii < len(cell_node_map.values)) if ii is not None else False for ii in ng_index] # Select only the points owned by this rank From e45542c7cafd6e21ccd730c2649c60dde91a58b1 Mon Sep 17 00:00:00 2001 From: Jack Betteridge Date: Thu, 3 Oct 2024 23:15:20 +0100 Subject: [PATCH 10/20] Remove even more parallel communication --- ngsPETSc/utils/firedrake/meshes.py | 38 ++++++++++++++++++++---------- 1 file changed, 26 insertions(+), 12 deletions(-) diff --git a/ngsPETSc/utils/firedrake/meshes.py b/ngsPETSc/utils/firedrake/meshes.py index 78801b2..eb4da5b 100644 --- a/ngsPETSc/utils/firedrake/meshes.py +++ b/ngsPETSc/utils/firedrake/meshes.py @@ -1,6 +1,6 @@ ''' This module contains all the functions related to wrapping NGSolve meshes to Firedrake -We adopt the same docstring conventiona as the Firedrake project, since this part of +We adopt the same docstring conventions as the Firedrake project, since this part of the package will only be used in combination with Firedrake. ''' try: @@ -164,26 +164,43 @@ def curveField(self, order, tol=1e-8): ) # Curve the mesh on rank 0 only - # JBTODO: (why?) if self.comm.rank == 0: + # Construct numpy arrays for physical domain data + physical_space_points = np.ndarray( + (ng_dimension, reference_space_points.shape[0], geom_dim) + ) + curved_space_points = np.ndarray( + (ng_dimension, reference_space_points.shape[0], geom_dim) + ) self.netgen_mesh.CalcElementMapping(reference_space_points, physical_space_points) self.netgen_mesh.Curve(order) self.netgen_mesh.CalcElementMapping(reference_space_points, curved_space_points) curved = ng_element().NumPy()["curved"] + # Broadcast a boolean array identifying curved cells + curved = self.comm.bcast(curved, root=0) + physical_space_points = physical_space_points[curved] + curved_space_points = curved_space_points[curved] else: - curved = np.array((ng_dimension, 1)) + curved = self.comm.bcast(None, root=0) + # Construct numpy arrays as buffers to receive physical domain data + ncurved = np.sum(curved) + physical_space_points = np.ndarray( + (ncurved, reference_space_points.shape[0], geom_dim) + ) + curved_space_points = np.ndarray( + (ncurved, reference_space_points.shape[0], geom_dim) + ) - # Broadcast curving data - physical_space_points = self.comm.bcast(physical_space_points, root=0) - curved_space_points = self.comm.bcast(curved_space_points, root=0) - curved = self.comm.bcast(curved, root=0) + # Broadcast curved cell point data + self.comm.Bcast(physical_space_points, root=0) + self.comm.Bcast(curved_space_points, root=0) cell_node_map = new_coordinates.cell_node_map() # Select only the points in curved cells - physical_space_points = physical_space_points[curved] - curved_space_points = curved_space_points[curved] barycentres = np.average(physical_space_points, axis=1) ng_index = [*map(lambda x: self.locate_cell(x, tolerance=0.01), barycentres)] + + # Select only the indices of points owned by this rank owned = [(0 <= ii < len(cell_node_map.values)) if ii is not None else False for ii in ng_index] # Select only the points owned by this rank @@ -196,11 +213,8 @@ def curveField(self, order, tol=1e-8): pyop2_index = [] for ngidx in ng_index: pyop2_index.extend(cell_node_map.values[ngidx]) - np.array(pyop2_index) # Find the correct coordinate permutation for each cell - # JBTODO: This should be moved to the next loop if we have to loop - # over cells any way to actually perform the permutation permutation = find_permutation( physical_space_points, new_coordinates.dat.data[pyop2_index].reshape(physical_space_points.shape), From 7b4f590a8ecaa235b8ccd8be8583d6f2af4b78f9 Mon Sep 17 00:00:00 2001 From: Jack Betteridge Date: Thu, 3 Oct 2024 23:39:30 +0100 Subject: [PATCH 11/20] Add PETSc event decorators --- ngsPETSc/utils/firedrake/meshes.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/ngsPETSc/utils/firedrake/meshes.py b/ngsPETSc/utils/firedrake/meshes.py index eb4da5b..5616fe6 100644 --- a/ngsPETSc/utils/firedrake/meshes.py +++ b/ngsPETSc/utils/firedrake/meshes.py @@ -86,9 +86,10 @@ def _slow_cdist(XA, XB): if not HAVE_SCIPY: - cdist = _slow_cdist + cdist = PETSc.Log.EventDecorator()(_slow_cdist) +@PETSc.Log.EventDecorator() def find_permutation(points_a, points_b, tol=1e-5): """ Find all permutations between a list of two sets of points. @@ -121,6 +122,7 @@ def find_permutation(points_a, points_b, tol=1e-5): return permutation +@PETSc.Log.EventDecorator() def curveField(self, order, tol=1e-8): ''' This method returns a curved mesh as a Firedrake function. From 3831f7c248266a781e32a37f9ce1db235ce9ffee Mon Sep 17 00:00:00 2001 From: Jack Betteridge Date: Fri, 4 Oct 2024 11:29:59 +0100 Subject: [PATCH 12/20] Prefer np.zeros over np.ndarray --- ngsPETSc/utils/firedrake/meshes.py | 16 ++++------------ 1 file changed, 4 insertions(+), 12 deletions(-) diff --git a/ngsPETSc/utils/firedrake/meshes.py b/ngsPETSc/utils/firedrake/meshes.py index 2939e40..7d492af 100644 --- a/ngsPETSc/utils/firedrake/meshes.py +++ b/ngsPETSc/utils/firedrake/meshes.py @@ -160,21 +160,13 @@ def curveField(self, order, tol=1e-8, CG=False): ref.append(pt) reference_space_points = np.array(ref) - # Map to the physical domain - physical_space_points = np.ndarray( - (ng_dimension, reference_space_points.shape[0], geom_dim) - ) - curved_space_points = np.ndarray( - (ng_dimension, reference_space_points.shape[0], geom_dim) - ) - # Curve the mesh on rank 0 only if self.comm.rank == 0: # Construct numpy arrays for physical domain data - physical_space_points = np.ndarray( + physical_space_points = np.zeros( (ng_dimension, reference_space_points.shape[0], geom_dim) ) - curved_space_points = np.ndarray( + curved_space_points = np.zeros( (ng_dimension, reference_space_points.shape[0], geom_dim) ) self.netgen_mesh.CalcElementMapping(reference_space_points, physical_space_points) @@ -189,10 +181,10 @@ def curveField(self, order, tol=1e-8, CG=False): curved = self.comm.bcast(None, root=0) # Construct numpy arrays as buffers to receive physical domain data ncurved = np.sum(curved) - physical_space_points = np.ndarray( + physical_space_points = np.zeros( (ncurved, reference_space_points.shape[0], geom_dim) ) - curved_space_points = np.ndarray( + curved_space_points = np.zeros( (ncurved, reference_space_points.shape[0], geom_dim) ) From bfa64872473732e4393e45af36c1c1a9466bd8a6 Mon Sep 17 00:00:00 2001 From: Umberto Zerbinati Date: Fri, 4 Oct 2024 12:08:26 +0100 Subject: [PATCH 13/20] Update install.rst to include @JDBetteridge as an author Thanks @JDBetteridge for solving parallel dead lock and speeding up the code! --- docs/src/install.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/install.rst b/docs/src/install.rst index d0a336a..cff476a 100644 --- a/docs/src/install.rst +++ b/docs/src/install.rst @@ -90,7 +90,7 @@ We are now finally ready to install ngsPETSc: Authors ---------- -Patrick E. Farrell, Stefano Zampini, Umberto Zerbinati +Jack Betteridge, Patrick E. Farrell, Stefano Zampini, Umberto Zerbinati License --------------- From d986ee0765bf7a9824e12e50322d79caaf5ddff4 Mon Sep 17 00:00:00 2001 From: Umberto Zerbinati Date: Fri, 4 Oct 2024 12:22:37 +0100 Subject: [PATCH 14/20] Different tollerance for permutation and barycenter location --- ngsPETSc/utils/firedrake/meshes.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/ngsPETSc/utils/firedrake/meshes.py b/ngsPETSc/utils/firedrake/meshes.py index 7d492af..71d80b4 100644 --- a/ngsPETSc/utils/firedrake/meshes.py +++ b/ngsPETSc/utils/firedrake/meshes.py @@ -123,11 +123,13 @@ def find_permutation(points_a, points_b, tol=1e-5): @PETSc.Log.EventDecorator() -def curveField(self, order, tol=1e-8, CG=False): +def curveField(self, order, perm_tol=1e-8, bary_tol=1e-1, CG=False): ''' This method returns a curved mesh as a Firedrake function. :arg order: the order of the curved mesh + :arg perm_tol: tolerance used to construct the permutation of the reference element. + :arg bary_tol: tolerance used to locate the cell a point belongs to. ''' # Check if the mesh is a surface mesh or two dimensional mesh @@ -195,7 +197,7 @@ def curveField(self, order, tol=1e-8, CG=False): # Select only the points in curved cells barycentres = np.average(physical_space_points, axis=1) - ng_index = [*map(lambda x: self.locate_cell(x, tolerance=0.01), barycentres)] + ng_index = [*map(lambda x: self.locate_cell(x, tolerance=bary_tol), barycentres)] # Select only the indices of points owned by this rank owned = [(0 <= ii < len(cell_node_map.values)) if ii is not None else False for ii in ng_index] @@ -215,7 +217,7 @@ def curveField(self, order, tol=1e-8, CG=False): permutation = find_permutation( physical_space_points, new_coordinates.dat.data[pyop2_index].reshape(physical_space_points.shape), - tol=tol + tol=perm_tol ) # Apply the permutation to each cell in turn From 0ecbfe7788aa76a8d22da4e9ad77335a6fc5bade Mon Sep 17 00:00:00 2001 From: Umberto Zerbinati Date: Fri, 4 Oct 2024 12:27:52 +0100 Subject: [PATCH 15/20] Scipy dependencies in setup.py --- setup.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/setup.py b/setup.py index fd9f332..4cc084f 100644 --- a/setup.py +++ b/setup.py @@ -2,7 +2,9 @@ with open("README.md", "r", encoding = "utf-8") as fh: long_description = fh.read() - + +scipy = '' if 'NGSPETSC_NO_SCIPY' in os.environ else 'scipy' + if 'NGSPETSC_NO_INSTALL_REQUIRED' in os.environ: install_requires = [] elif 'NGS_FROM_SOURCE' in os.environ: @@ -10,6 +12,7 @@ 'petsc4py', 'mpi4py', 'numpy', + scipy, 'pytest', #For testing 'pylint', #For formatting ] @@ -19,6 +22,7 @@ 'ngsolve', 'petsc4py', 'mpi4py', + scipy, 'pytest', #For testing 'pylint', #For formatting ] From 9824948d7560563af172015ef6fa5adf6149ae59 Mon Sep 17 00:00:00 2001 From: Umberto Zerbinati Date: Fri, 4 Oct 2024 12:31:51 +0100 Subject: [PATCH 16/20] Update meshes.py --- ngsPETSc/utils/firedrake/meshes.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/ngsPETSc/utils/firedrake/meshes.py b/ngsPETSc/utils/firedrake/meshes.py index 71d80b4..cdfe6fb 100644 --- a/ngsPETSc/utils/firedrake/meshes.py +++ b/ngsPETSc/utils/firedrake/meshes.py @@ -123,13 +123,15 @@ def find_permutation(points_a, points_b, tol=1e-5): @PETSc.Log.EventDecorator() -def curveField(self, order, perm_tol=1e-8, bary_tol=1e-1, CG=False): +def curveField(self, order, perm_tol=1e-8, bary_tol=1e-1, cg_field=False): ''' This method returns a curved mesh as a Firedrake function. - :arg order: the order of the curved mesh + :arg order: the order of the curved mesh. :arg perm_tol: tolerance used to construct the permutation of the reference element. :arg bary_tol: tolerance used to locate the cell a point belongs to. + :arg cg_field: return a CG function field representing the mesh, rather than the + default DG field. ''' # Check if the mesh is a surface mesh or two dimensional mesh @@ -141,7 +143,7 @@ def curveField(self, order, perm_tol=1e-8, bary_tol=1e-1, CG=False): geom_dim = self.geometric_dimension() # Construct the mesh as a Firedrake function - if CG: + if cg_field: firedrake_space = fd.VectorFunctionSpace(self, "CG", order) else: low_order_element = self.coordinates.function_space().ufl_element().sub_elements[0] From 14f27190e1b13e03184024abd1e174cdf386e719 Mon Sep 17 00:00:00 2001 From: Jack Betteridge Date: Fri, 4 Oct 2024 14:16:47 +0100 Subject: [PATCH 17/20] Move stuff to pyproject.toml --- pyproject.toml | 38 ++++++++++++++++++++++++++++++++++++++ setup.py | 30 ++++-------------------------- 2 files changed, 42 insertions(+), 26 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index f1bfca4..269d8c6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,3 +1,41 @@ +[project] +name = "ngsPETSc" +version = "0.0.5" +description = "NGSolve/Netgen interface to PETSc." +readme = "README.md" +authors = [ + {name = "Umberto Zerbinati", email = "umberto.zerbinati@maths.ox.ac.uk"}, + {name = "Patrick E. Farrell", email = "patrick.farrell@maths.ox.ac.uk"}, + {name = "Stefano Zampini", email = "stefano.zampini@kaust.edu.sa"}, + {name = "Jack Betteridge", email = "J.Betteridge@imperial.ac.uk"}, +] +maintainers = [ + {name = "Umberto Zerbinati", email = "umberto.zerbinati@maths.ox.ac.uk"}, +] +license = {file = "LICENSE.txt"} +dependencies = [ + "mpi4py", + "numpy", + "scipy", +] +requires-python = ">=3.8" +classifiers = [ + "Programming Language :: Python :: 3", + "License :: OSI Approved :: MIT License", + "Operating System :: OS Independent", + "Development Status :: 3 - Alpha", +] + +[project.urls] +Documentation = "https://ngspetsc.readthedocs.io/en/latest/" +Repository = "https://github.com/NGSolve/ngsPETSc" + +[project.optional-dependencies] +dev = [ + "pytest", + "pylint", +] + [build-system] requires = ['setuptools>=42'] build-backend = 'setuptools.build_meta' diff --git a/setup.py b/setup.py index 4cc084f..5ef6ff9 100644 --- a/setup.py +++ b/setup.py @@ -1,10 +1,5 @@ import setuptools, os -with open("README.md", "r", encoding = "utf-8") as fh: - long_description = fh.read() - -scipy = '' if 'NGSPETSC_NO_SCIPY' in os.environ else 'scipy' - if 'NGSPETSC_NO_INSTALL_REQUIRED' in os.environ: install_requires = [] elif 'NGS_FROM_SOURCE' in os.environ: @@ -12,7 +7,7 @@ 'petsc4py', 'mpi4py', 'numpy', - scipy, + 'scipy', 'pytest', #For testing 'pylint', #For formatting ] @@ -22,29 +17,12 @@ 'ngsolve', 'petsc4py', 'mpi4py', - scipy, + 'scipy', 'pytest', #For testing 'pylint', #For formatting ] setuptools.setup( - name = "ngsPETSc", - version = "0.0.5", - author = "Umberto Zerbinati", - author_email = "umberto.zerbinati@maths.ox.ac.uk", - description = "NGSolve/Netgen interface to PETSc.", - long_description = long_description, - long_description_content_type = "text/markdown", - url = "", - project_urls = { - }, - classifiers = [ - "Programming Language :: Python :: 3", - "License :: OSI Approved :: MIT License", - "Operating System :: OS Independent", - ], - packages=["ngsPETSc", "ngsPETSc.utils", "ngsPETSc.utils.firedrake", "ngsPETSc.utils.ngs"], - python_requires = ">=3.8", - install_requires=install_requires - + install_requires=install_requires, + packages=["ngsPETSc", "ngsPETSc.utils", "ngsPETSc.utils.firedrake", "ngsPETSc.utils.ngs"] ) From ecb0a2d60461dd3790e2f30ab167b815a02e15e7 Mon Sep 17 00:00:00 2001 From: Jack Betteridge Date: Fri, 4 Oct 2024 14:20:42 +0100 Subject: [PATCH 18/20] Explicit names --- ngsPETSc/utils/firedrake/meshes.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/ngsPETSc/utils/firedrake/meshes.py b/ngsPETSc/utils/firedrake/meshes.py index cdfe6fb..c6cb755 100644 --- a/ngsPETSc/utils/firedrake/meshes.py +++ b/ngsPETSc/utils/firedrake/meshes.py @@ -123,14 +123,14 @@ def find_permutation(points_a, points_b, tol=1e-5): @PETSc.Log.EventDecorator() -def curveField(self, order, perm_tol=1e-8, bary_tol=1e-1, cg_field=False): +def curveField(self, order, permutation_tol=1e-8, location_tol=1e-1, cg_field=False): ''' This method returns a curved mesh as a Firedrake function. :arg order: the order of the curved mesh. - :arg perm_tol: tolerance used to construct the permutation of the reference element. - :arg bary_tol: tolerance used to locate the cell a point belongs to. - :arg cg_field: return a CG function field representing the mesh, rather than the + :arg permutation_tol: tolerance used to construct the permutation of the reference element. + :arg location_tol: tolerance used to locate the cell a point belongs to. + :arg cg_field: return a CG function field representing the mesh, rather than the default DG field. ''' @@ -199,7 +199,7 @@ def curveField(self, order, perm_tol=1e-8, bary_tol=1e-1, cg_field=False): # Select only the points in curved cells barycentres = np.average(physical_space_points, axis=1) - ng_index = [*map(lambda x: self.locate_cell(x, tolerance=bary_tol), barycentres)] + ng_index = [*map(lambda x: self.locate_cell(x, tolerance=location_tol), barycentres)] # Select only the indices of points owned by this rank owned = [(0 <= ii < len(cell_node_map.values)) if ii is not None else False for ii in ng_index] @@ -219,7 +219,7 @@ def curveField(self, order, perm_tol=1e-8, bary_tol=1e-1, cg_field=False): permutation = find_permutation( physical_space_points, new_coordinates.dat.data[pyop2_index].reshape(physical_space_points.shape), - tol=perm_tol + tol=permutation_tol ) # Apply the permutation to each cell in turn From 12b5ddd2470c91988a6bba47fa90f8ab57edf0ae Mon Sep 17 00:00:00 2001 From: Jack Betteridge Date: Fri, 4 Oct 2024 14:25:56 +0100 Subject: [PATCH 19/20] Update hierarchies for interface change --- ngsPETSc/utils/firedrake/hierarchies.py | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/ngsPETSc/utils/firedrake/hierarchies.py b/ngsPETSc/utils/firedrake/hierarchies.py index 286c88d..2aaae93 100644 --- a/ngsPETSc/utils/firedrake/hierarchies.py +++ b/ngsPETSc/utils/firedrake/hierarchies.py @@ -1,5 +1,5 @@ ''' -This module contains all the functions related +This module contains all the functions related ''' try: import firedrake as fd @@ -204,7 +204,7 @@ def NetgenHierarchy(mesh, levs, flags): order = flagsUtils(flags, "degree", 1) if isinstance(order, int): order= [order]*(levs+1) - tol = flagsUtils(flags, "tol", 1e-8) + permutation_tol = flagsUtils(flags, "tol", 1e-8) refType = flagsUtils(flags, "refinement_type", "uniform") optMoves = flagsUtils(flags, "optimisation_moves", False) snap = flagsUtils(flags, "snap_to", "geometry") @@ -221,7 +221,11 @@ def NetgenHierarchy(mesh, levs, flags): raise RuntimeError("Cannot refine parallel overlapped meshes ") #We curve the mesh if order[0]>1: - ho_field = mesh.curve_field(order=order[0], tol=tol, CG=cg) + ho_field = mesh.curve_field( + order=order[0], + permutation_tol=permutation_tol, + cg_field=cg + ) mesh = fd.Mesh(ho_field,distribution_parameters=params, comm=comm) meshes += [mesh] cdm = meshes[-1].topology_dm @@ -248,8 +252,11 @@ def NetgenHierarchy(mesh, levs, flags): #We curve the mesh if order[l+1] > 1: if snap == "geometry": - mesh = fd.Mesh(mesh.curve_field(order=order[l+1], tol=tol), - distribution_parameters=params, comm=comm) + mesh = fd.Mesh( + mesh.curve_field(order=order[l+1], permutation_tol=permutation_tol), + distribution_parameters=params, + comm=comm + ) elif snap == "coarse": mesh = snapToCoarse(ho_field, mesh, order[l+1], snap_smoothing, cg) meshes += [mesh] From 4874b86c398b10a6178735d6e0e482d74cce646c Mon Sep 17 00:00:00 2001 From: Jack Betteridge Date: Fri, 4 Oct 2024 14:30:30 +0100 Subject: [PATCH 20/20] Remove slow cdist as scipy is now a hard dependency --- ngsPETSc/utils/firedrake/meshes.py | 19 +------------------ 1 file changed, 1 insertion(+), 18 deletions(-) diff --git a/ngsPETSc/utils/firedrake/meshes.py b/ngsPETSc/utils/firedrake/meshes.py index c6cb755..4fc10a9 100644 --- a/ngsPETSc/utils/firedrake/meshes.py +++ b/ngsPETSc/utils/firedrake/meshes.py @@ -11,12 +11,7 @@ import numpy as np from petsc4py import PETSc - -try: - from scipy.spatial.distance import cdist - HAVE_SCIPY = True -except ImportError: - HAVE_SCIPY = False +from scipy.spatial.distance import cdist import netgen import netgen.meshing as ngm @@ -77,18 +72,6 @@ def refineMarkedElements(self, mark): raise NotImplementedError("No implementation for dimension other than 2 and 3.") -def _slow_cdist(XA, XB): - dist = np.zeros([len(XA), len(XB)]) - for ii, a in enumerate(XA): - for jj, b in enumerate(XB): - dist[ii, jj] = np.linalg.norm(b - a) - return dist - - -if not HAVE_SCIPY: - cdist = PETSc.Log.EventDecorator()(_slow_cdist) - - @PETSc.Log.EventDecorator() def find_permutation(points_a, points_b, tol=1e-5): """ Find all permutations between a list of two sets of points.