Added new mesh creation framework

cardmech/conf/elasticity.conf

@@ -45,8 +45,8 @@ PressureIterations = 1
 
 DebugPressure = 1
 
-MechLevel = 0
-MechPolynomialDegree = 1
+MechLevel = 1
+MechPolynomialDegree = 2
 
 WithPrestress=false
 PrestressSteps=5
@@ -59,6 +59,7 @@ PrestressSteps=5
 #       Meshes       #
 ######################
 #Mesh = BenchmarkBeam2DTet
+Mesh = QuarterEllipsoid
 Mesh = OrientedEllipsoid
 #Mesh = KovachevaHeart
 #Mesh = VolumeBlock3DTet
@@ -178,6 +179,12 @@ StVenantMat_Mu = 8000
 # Guccione Material
 #------------------------
 
+# Benchmark Parameter
+GuccioneMat_C=10000
+GuccioneMat_bf=1
+GuccioneMat_bs=1
+GuccioneMat_bfs=1
+
 # Benchmark Parameter
 GuccioneMat_C=2000
 GuccioneMat_bf=80
@@ -235,7 +242,7 @@ HolzapfelMat_bfs = 9.466
 # Pressure Parameters#
 ######################
 #ConstantPressure = -10
-LVPressure = 15000
+LVPressure = 10000
 RVPressure = 0
 RAPressure = 0
 LAPressure = 0

tools/geometries/utility/geo.py

@@ -64,13 +64,13 @@ def write_properties_oriented(orientations, activation, transform=None):
 
 
 class Geometry:
+
     def __init__(self, name, path="", withdata=False):
         self.name = name
         self.path = path
         self.points = {}
         self.subdomains = {}
         self.cells = {}
-        self.bnd = {}
         self.faces = {}
         self.vdata = {}
         self.cdata = {}
@@ -108,6 +108,7 @@ class Geometry:
 
     def remove_face(self, face_index):
         del self.faces[face_index]
+        del self.bnd[face_index]
 
     def init_data(self):
         self.vdata = {i: [] for i in self.points.keys()}

tools/geometries/utility/mesh.py

+import itertools
+
+precision = 6
+
+default_vdata = [-1, 0, 0]
+default_cdata = [1, 0, 0, 0, 1, 0, 0, 0, 1]
+
+
+class cell:
+    def __init__(self, i, point_indices, face_size):
+        self.index = i
+        self._indices = point_indices
+        self._fsize = face_size
+        self.subdomain = 0
+
+    def __len__(self):
+        return len(self._indices)
+
+    def faces(self):
+        return itertools.combinations(self._indices, self._fsize)
+
+    def corners(self):
+        return self._indices
+
+
+class face:
+    def __init__(self, lcell_index, point_indices):
+        self.left = lcell_index
+        self.right = -1
+        self.bnd = -1
+        self._indices = point_indices
+
+    def set_rside(self, rcell_index):
+        self.right = rcell_index
+
+    def on_bnd(self):
+        return self.right < 0
+
+    def __len__(self):
+        return len(self._indices)
+
+    def corners(self):
+        return self._indices
+
+
+def matches(meshobject, points):
+    return all(p in points for p in meshobject.corners())
+
+
+class Mesh:
+    def __init__(self, name, path="", withdata=False):
+        self.name = name
+        self.path = path
+        self.points = {}
+        self.cells = {}
+        self.faces = {}
+        self.vdata = {}
+        self.cdata = {}
+        self.is_data_geo = withdata
+
+    def add_point(self, point):
+        index = len(self.points)
+        self.points[index] = [round(p, precision) for p in point]
+        return index
+
+    def add_points(self, points):
+        for point in points:
+            self.add_point(point)
+
+    def add_cell(self, p_indices, fsize):
+        index = len(self.cells)
+        c = cell(index, p_indices, fsize)
+        self.cells[index] = c
+
+        for f in c.faces():
+            self.add_face(index, f)
+
+        return index
+
+    def add_cells(self, cells, fsize):
+        for c in cells:
+            self.add_cell(c, fsize)
+
+    def remove_cell(self, cell_index):
+        del self.cells[cell_index]
+
+    def add_face(self, lcell_index, p_indices):
+        existing_f = [f for f in self.faces.values() if matches(f, p_indices)]
+        if len(existing_f) > 0:
+            existing_f[0].right = lcell_index
+        else:
+            index = len(self.faces)
+            f = face(lcell_index, p_indices)
+            self.faces[index] = f
+
+    def remove_face(self, face_index):
+        del self.faces[face_index]
+
+    def init_data(self):
+        self.vdata = {i: [] for i in self.points.keys()}
+        self.cdata = {i: [] for i in self.cells.keys()}
+
+    def add_vdata(self, p_index, data):
+        if p_index not in self.vdata:
+            self.vdata[p_index] = []
+        for d in data:
+            self.vdata[p_index].append(d)
+
+    def add_vdata_array(self, v_data):
+        for i in range(len(v_data)):
+            self.add_vdata(i, v_data[i])
+
+    def fill_vdata(self, data):
+        for i, vd in self.vdata.items():
+            if len(vd) == 0:
+                self.add_vdata(i, data)
+
+    def add_cdata(self, c_index, data):
+        if c_index not in self.cdata:
+            self.cdata[c_index] = []
+        for d in data:
+            self.cdata[c_index].append(d)
+
+    def add_cdata_array(self, c_data):
+        for i in range(len(c_data)):
+            self.add_cdata(i, c_data[i])
+
+    def fill_cdata(self, data):
+        for i, cd in self.cdata.items():
+            if len(cd) == 0:
+                self.add_cdata(i, data)
+
+    def add_subdomain(self, c_index, sd):
+        self.cells[c_index].sd = sd
+
+    def add_bnd_condition(self, f_index, bnd):
+        self.faces[f_index].bnd = bnd
+
+    def fill_data(self, data="12 1 0 0 0 1 0 0 0 1 -10 0 0"):
+        for i in self.points.keys():
+            self.vdata[i] = data
+        for i in self.cells.keys():
+            self.cdata[i] = data
+
+    def point_count(self):
+        return len(self.points)
+
+    def cell_count(self):
+        return len(self.cells)
+
+    def face_count(self):
+        return len(self.faces)
+
+    def get_points(self):
+        return [p for p in self.points.values()]
+
+    def get_cells(self):
+        return [c for c in self.cells.values()]
+
+    def bnd_faces(self):
+        return [f for f in self.faces.values() if f.on_bnd()]
+
+    def corners(self, meshobject):
+        return [self.points[i] for i in meshobject.corners()]
+
+    def consistency(self):
+        missing__indices = [i for i in self.points if i not in self.vdata]
+        for i in missing__indices:
+            self.add_vdata(i, default_vdata)
+
+        missing__indices = [i for i in self.cells if i not in self.cdata]
+        for i in missing__indices:
+            self.add_cdata(i, default_cdata)
+
+    def save(self):
+        self.consistency()
+
+        with open(str(self.path + self.name + '.geo'), 'w') as output_file:
+            output_file.write("POINTS:\n")
+            output_file.write(self.print_points())
+
+            output_file.write("CELLS:\n")
+            output_file.write(self.print_cells())
+
+            output_file.write("FACES:\n")
+            output_file.write(self.print_faces())
+
+            if self.is_data_geo:
+                output_file.write("VDATA:\n")
+                output_file.write(self.print_pointdata())
+
+                output_file.write("CDATA:\n")
+                output_file.write(self.print_celldata())
+
+    def print_points(self):
+        points_as_string = ""
+        for point in self.points.values():
+            for p in point:
+                points_as_string += str(p) + " "
+            points_as_string = points_as_string[:-1] + "\n"
+        return points_as_string
+
+    def print_cells(self):
+        cells_as_string = ""
+        for i, cell in self.cells.items():
+            cells_as_string += str(len(cell)) + " " + str(cell.subdomain)
+            for p in cell.corners():
+                cells_as_string += " " + str(p)
+            cells_as_string += "\n"
+        return cells_as_string
+
+    def print_faces(self):
+        faces_as_string = ""
+        for i, face in self.faces.items():
+            if not face.on_bnd():
+                continue
+            faces_as_string += str(len(face)) + " " + str(max(face.bnd, 0))
+            for p in face.corners():
+                faces_as_string += " " + str(p)
+            faces_as_string += "\n"
+        return faces_as_string
+
+    def print_pointdata(self):
+        data_as_string = ""
+        for data in self.vdata.values():
+            data_as_string += str(len(data))
+            for d in data:
+                data_as_string += " " + str(d)
+            data_as_string += "\n"
+        return data_as_string
+
+    def print_celldata(self):
+        data_as_string = ""
+        for data in self.cdata.values():
+            data_as_string += str(len(data))
+            for d in data:
+                data_as_string += " " + str(d)
+            data_as_string += "\n"
+        return data_as_string
+
+    def cell_midpoint(self, cell_id):
+        cell = self.cells[cell_id]
+        midpoint = [0 for _ in range(len(self.points[0]))]
+        dim = len(midpoint)
+        for p_id in cell:
+            p = self.points[p_id]
+            for i in range(dim):
+                midpoint[i] += p[i] / len(cell)
+
+        return midpoint
+
+    def find_face(self, point_indices):
+        for face in self.faces:
+            if matches(face, point_indices):
+                return face
+        raise KeyError("Face not found")

tools/geometries/vtu_to_geo.py

@@ -382,6 +382,7 @@ def convert_2LayerEllipsoid():
     activationByPointList(geometry,True)
     finish(geometry)
 
+
 def convert_ellipsoid():
     vtu = VtuGeometry("vtk/", "ellipsoid", "OrientedEllipsoid")
     vtu.rescale(1000.0)
@@ -399,28 +400,27 @@ def convert_ellipsoid():
     for fi in faces_to_remove:
         geometry.remove_face(fi)
     activation_ventricle(geometry)
+    print("Cells: ", len(geometry.cells))
     finish(geometry)
 
-    #vtu = VtuGeometry("vtk/", "ellipsoid", "EllipsoidExcitationApex")
-    #vtu.rescale(1000.0)
-    #mat_map = vtu.material_map
-    #for k in range(4):
-        #mat_map[30 + k] = 0
-       # mat_map[10 + k] = 230
-   # vtu.remap(mat_map)
-
-    #geometry = vtu.convert(1, 1, [], [2, 3, 4])
-    #faces_to_remove = []
-    #for fi in geometry.faces:
-       # if geometry.bnd[fi] < 100:
-           # faces_to_remove.append(fi)
-    #for fi in faces_to_remove:
-       # geometry.remove_face(fi)
-    #no_activaton(geometry)
-    #finish(geometry)
-    
-    
-    
+    # vtu = VtuGeometry("vtk/", "ellipsoid", "EllipsoidExcitationApex")
+    # vtu.rescale(1000.0)
+    # mat_map = vtu.material_map
+    # for k in range(4):
+    # mat_map[30 + k] = 0
+    # mat_map[10 + k] = 230
+    # vtu.remap(mat_map)
+
+    # geometry = vtu.convert(1, 1, [], [2, 3, 4])
+    # faces_to_remove = []
+    # for fi in geometry.faces:
+    # if geometry.bnd[fi] < 100:
+    # faces_to_remove.append(fi)
+    # for fi in faces_to_remove:
+    # geometry.remove_face(fi)
+    # no_activaton(geometry)
+    # finish(geometry)
+
     vtu = VtuGeometry("vtk/", "ellipsoid", "PlainEllipsoid")
     vtu.rescale(1000.0)
     mat_map = vtu.material_map
@@ -439,13 +439,50 @@ def convert_ellipsoid():
     activation_ventricle(geometry)
     finish(geometry)
 
+
+def quarter_ellipsoid_bnd(geometry):
+    # geometry.bnd.clear()
+    # geometry.faces.clear()
+
+    for ci in geometry.cells:
+        for face in geometry.cell_faces(ci):
+            if all(geometry.points[pi][0] == 0.0 for pi in face):
+                fi = geometry.add_face(face)
+                geometry.add_bnd_condition(fi, 100)
+            elif all(geometry.points[pi][1] == 0.0 for pi in face):
+                fi = geometry.add_face(face)
+                geometry.add_bnd_condition(fi, 101)
+
+            # elif all(geometry.points[pi][2] == 5.0 for pi in face):
+            #    fi = geometry.add_face(face)
+            #    geometry.add_bnd_condition(fi, 199)
+
+
+def convert_quarter_ellipsoid():
+    vtu = VtuGeometry("vtk/", "QuarterEllipsoid", "QuarterEllipsoid")
+    vtu.rescale(1.0)
+    mat_map = vtu.material_map
+    for k in range(4):
+        mat_map[30 + k] = 0
+        mat_map[10 + k] = 230
+    vtu.remap(mat_map)
+
+    geometry = vtu.convert(0, -1, [], [1, 2, 3])
+    quarter_ellipsoid_bnd(geometry)
+
+    activation_ventricle(geometry)
+    print("Cells: ", len(geometry.cells))
+    finish(geometry)
+
+
 def robin_bnd_ellispoid(vtuconv, geometry, surface_data):
     surface_points = []
     surface_indices = []
     with open(surface_data) as datafile:
         lines = [line.split(',') for line in datafile.readlines()]
         for line in lines[1:]:
-            surface_points.append([vtuconv.scale * float(line[-3]),  vtuconv.scale * float(line[-2]),  vtuconv.scale * float(line[-1])])
+            surface_points.append(
+                [vtuconv.scale * float(line[-3]), vtuconv.scale * float(line[-2]), vtuconv.scale * float(line[-1])])
 
     spoints = len(surface_points)
     for _ in range(spoints):
@@ -490,12 +527,13 @@ def transform_robin_ellipsoid():
 
 if __name__ == "__main__":
     # convert_T4()
-    convert_2LayerEllipsoid()
-    #convert_plainEllipsoid()
-    #convert_ellipsoid()
-    #convert_ventricle()
-    #convert_heart_kovacheva()
-
-    #transform_robin_ellipsoid()
+    # convert_2LayerEllipsoid()
+    # convert_plainEllipsoid()
+    convert_ellipsoid()
+    convert_quarter_ellipsoid()
+    # convert_ventricle()
+    # convert_heart_kovacheva()
+
+    # transform_robin_ellipsoid()
     # convert_vtu("../../../IBT/EP_EK/EP_EK/source/", "Heart", "HeartReg", 0, -1, [], [9, 10, 11])
     # convert_vtu("../../../IBT/Geomtries/", "heartT4", "HeartT4", 1, -1, [], [2, 3, 4])

tools/geometries/vtu_to_mesh.py

+from vtkmodules.util.numpy_support import vtk_to_numpy
+from vtkmodules.vtkIOXML import vtkXMLUnstructuredGridReader
+
+from utility.mesh import Mesh
+
+
+class VtuMesh:
+    def __init__(self, vtu_file, geo_name, scale=1.0):
+        self.geo = Mesh(geo_name, "../../geo/", True)
+        self.bnd_data = []
+        self.scale = scale
+
+        # Read File
+        reader = vtkXMLUnstructuredGridReader()
+        reader.SetFileName(vtu_file + ".vtu")
+        reader.Update()
+        data = reader.GetOutput()
+
+        # Extract Points
+        nodes = vtk_to_numpy(data.GetPoints().GetData())
+        self.add_points(nodes)
+        # Extract Cells
+        self.add_cells(vtk_to_numpy(data.GetCells().GetData()))
+
+        self.point_data = data.GetPointData()
+        self.cell_data = data.GetCellData()
+
+    def add_points(self, nodes):
+        for node in nodes:
+            self.geo.add_point([p * self.scale for p in node])
+
+    def add_cells(self, cells):
+        t = cells.size
+        j, k = 0, 0
+        while j < t:
+            i = cells[j]
+            if i == 10 or i == 4:
+                self.geo.add_cell([p for p in cells[j + 1:j + 5]], 3)
+            elif i == 3:
+                self.bnd_data.append([p for p in cells[j + 1:j + i + 1]])
+            j += i + 1
+
+    def set_boundary_data(self, bnd_function):
+        for bnd_face in self.geo.bnd_faces():
+            indices = bnd_face.corners()
+            points = self.geo.corners(bnd_face)
+            bnd_face.bnd = bnd_function({indices[i]: points[i] for i in range(len(indices))})
+
+    def read_celldata(self, cell_indices):
+        for ci in cell_indices:
+            data = vtk_to_numpy(self.cell_data.GetAbstractArray(ci))
+            self.geo.add_cdata_array(data[:self.geo.cell_count()])
+
+    def read_boundary_data(self, bnd_index, bnd_map=lambda x: x):
+        bnd_values = vtk_to_numpy(self.cell_data.GetAbstractArray(bnd_index))
+
+        shift = self.geo.cell_count()
+        for fi in range(len(self.bnd_data)):
+            bnd_cond = bnd_map(bnd_values[shift + fi])
+            bnd_face = self.geo.find_face(self.bnd_data[fi])
+            bnd_face.bnd = bnd_cond
+
+    def save(self):
+        self.geo.save()
+
+
+quarter_ellipsoid_epicard = (
+    298, 299, 300, 301, 302, 303, 304, 305, 306, 307, 308, 309, 310, 311, 312, 313, 314, 315, 316, 317, 318, 319, 320,
+    321,
+    322, 323, 324, 325, 326, 327, 328, 329, 330, 331, 332, 333, 334, 335, 336, 337, 338, 339, 340, 341, 342, 343, 344,
+    345,
+    346, 347, 348, 349, 350, 351, 352, 353, 354, 355, 356, 357, 358, 359, 360, 361, 362, 363, 364, 365, 366, 367, 368,
+    369,
+    370, 371, 372, 373, 374, 375, 376, 377, 378, 379, 380, 381, 382, 383, 384, 385, 386, 387, 388, 389, 390, 391, 392,
+    393,
+    394, 395, 396, 397, 398, 399, 400, 401, 402, 403, 404, 405, 406, 407, 408, 409, 410, 411, 412, 413, 414, 415, 416,
+    417,
+    418, 419, 420, 421, 422, 423, 424, 425, 426, 427, 428, 429, 430, 431, 432, 433, 434, 435, 436, 437, 438, 439, 440,
+    441,
+    442, 443, 444, 445, 446, 447, 448, 449, 450, 451, 452, 453, 454, 455, 456, 457, 458, 459, 460, 461, 462, 463, 464,
+    465,
+    466, 467, 468, 469, 470, 471, 472, 473, 474, 475, 476, 477, 478, 479, 480, 481, 482, 483, 484, 485, 486, 487, 488,
+    489,
+    490, 491, 492, 493, 494, 495, 496, 497, 498, 499, 500, 501, 502, 503)
+
+
+def quarter_ellipsoid_bnd(points):
+    if all([p[0] == 0 for p in points.values()]):
+        return 100
+    if all([p[1] == 0 for p in points.values()]):
+        return 101
+    if all([p[2] == 5 for p in points.values()]):
+        return 199
+
+    if all([p in quarter_ellipsoid_epicard for p in points.keys()]):
+        return 0
+
+    return 230
+
+
+if __name__ == '__main__':
+    geo = VtuMesh("vtk/QuarterEllipsoid", "QuarterEllipsoid", True)
+    geo.set_boundary_data(quarter_ellipsoid_bnd)
+    geo.read_celldata([1, 2, 3])
+    geo.save()
+    print("Saved new geometry")