Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
59 changes: 0 additions & 59 deletions src/pulse/geometry.py
Original file line number Diff line number Diff line change
Expand Up @@ -158,65 +158,6 @@ def comm(self) -> MPI.Intracomm:
return self.mesh.comm


class BaseData(NamedTuple):
centroid: npt.NDArray[np.float64]
vector: npt.NDArray[np.float64]
normal: npt.NDArray[np.float64]


def compute_base_data(
mesh: dolfinx.mesh.Mesh,
facet_tags: dolfinx.mesh.MeshTags,
marker,
) -> BaseData:
"""Compute the centroid, vector and normal of the base

Parameters
----------
mesh : dolfinx.mesh.Mesh
The mesh
facet_tags : dolfinx.mesh.MeshTags
The facet tags
marker : _type_
Marker for the base

Returns
-------
BaseData
NamedTuple containing the centroid, vector and normal of the base
"""
base_facets = facet_tags.find(marker)
base_midpoints = mesh.comm.gather(
dolfinx.mesh.compute_midpoints(mesh, 2, base_facets),
root=0,
)
base_vector = np.zeros(3)
base_centroid = np.zeros(3)
base_normal = np.zeros(3)
if mesh.comm.rank == 0:
bm = np.concatenate(base_midpoints)
base_centroid = bm.mean(axis=0)
# print("Base centroid", len(base_midpoints))
base_points_centered = bm - base_centroid
u, s, vh = np.linalg.svd(base_points_centered)
base_normal = vh[-1, :]
# Initialize vector to be used for cross product
vector_init = np.array([0, 1, 0])

# If the normal is parallel to the initial vector, change the initial vector
if np.allclose(np.abs(base_normal), np.abs(vector_init)):
vector_init = np.array([0, 0, 1])

# Find two vectors in the plane, orthogonal to the normal
vector = np.cross(base_normal, vector_init)
base_vector = np.cross(base_normal, vector)

base_centroid = mesh.comm.bcast(base_centroid, root=0)
base_vector = mesh.comm.bcast(base_vector, root=0)
base_normal = mesh.comm.bcast(base_normal, root=0)
return BaseData(centroid=base_centroid, vector=base_vector, normal=base_normal)


# Note slots doesn't work due to https://github.com/python/cpython/issues/90562
@dataclass(kw_only=True)
class HeartGeometry(Geometry):
Expand Down
2 changes: 0 additions & 2 deletions tests/test_geometry.py
Original file line number Diff line number Diff line change
@@ -1,6 +1,5 @@
import math

import gmsh
import numpy as np
import pytest
import ufl
Expand Down Expand Up @@ -112,7 +111,6 @@ def test_HeartGeometry_lv(tmp_path):
# assert np.isclose(geo2.volume("ENDO"), endo_volume, atol=1e-7)


@pytest.mark.skipif(gmsh.__version__ == "4.14.0", reason="GMSH 4.14.0 has a bug with fuse")
def test_HeartGeometry_biv(tmp_path):
geo1 = cardiac_geometries.mesh.biv_ellipsoid(
outdir=tmp_path,
Expand Down