diff --git a/CHANGELOG.md b/CHANGELOG.md index 761b90c..2e2cad2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +## [0.10.0] - 2022-10-12 +## Added +- possibility to have centre of mass as division point for the octree + ## [0.9.0] - 2022-10-09 ## Added - mesh2d as for diagnostic purposes diff --git a/dsmc/octree.py b/dsmc/octree.py index dd40ed0..54a8062 100644 --- a/dsmc/octree.py +++ b/dsmc/octree.py @@ -1,6 +1,7 @@ import numpy as np import numpy.typing as npt from numba import njit +from enum import Enum from . import common as com fmin = np.finfo(float).min @@ -121,33 +122,32 @@ def _create_boxes(box): return [child_geo1, child_geo2, child_geo3, child_geo4, child_geo5, child_geo6, child_geo7, child_geo8] @njit -def _get_min_aspect_ratio(box, axis): - half = np.array([0.5*(box[i][1] - box[i][0]) for i in range(3)]) +def _get_min_aspect_ratio(box, axis, half): + half_loc = np.array([0.5*(half[i] - box[i][0]) for i in range(3)]) match axis: case 0: - return min(half[0] / half[1], half[0] / half[2]); + return min(half_loc[0] / half_loc[1], half_loc[0] / half_loc[2]); case 1: - return min(half[1] / half[0], half[1] / half[2]); + return min(half_loc[1] / half_loc[0], half_loc[1] / half_loc[2]); case 2: - return min(half[2] / half[1], half[2] / half[0]); + return min(half_loc[2] / half_loc[1], half_loc[2] / half_loc[0]); @njit -def _devide(box, axis): - half = 0.5*(box[axis][0] + box[axis][1]) +def _devide(box, axis, half): box1 = np.copy(box) box2 = np.copy(box) box1[axis][0] = box[axis][0] - box1[axis][1] = half + box1[axis][1] = half[axis] - box2[axis][0] = half + box2[axis][0] = half[axis] box2[axis][1] = box[axis][1] return (box1, box2) @njit -def _create_combined_boxes(box, min_aspect_ratio): +def _create_combined_boxes(box, min_aspect_ratio, half): boxes = np.empty((15, 3, 2)) boxes[0] = box N = 0 @@ -155,9 +155,9 @@ def _create_combined_boxes(box, min_aspect_ratio): q = 1 for i in range(3): - if _get_min_aspect_ratio(box, i) > min_aspect_ratio: + if _get_min_aspect_ratio(box, i, half) > min_aspect_ratio: for b in range(Nold, Nold + 2**N): - new_boxes = _devide(boxes[b], i) + new_boxes = _devide(boxes[b], i, half) boxes[q] = new_boxes[0] boxes[q + 1] = new_boxes[1] q += 2 @@ -171,6 +171,20 @@ def _create_combined_boxes(box, min_aspect_ratio): new_boxes[b] = boxes[Nold + b] return new_boxes + +@njit +def _get_centre_of_mass(permutations, positions, offset, n_elements): + com = np.zeros(3) + + for i in range(offset, offset + n_elements): + p = permutations[i] + com += positions[p] + + return com / float(n_elements) + +class Type(Enum): + COV = 0 + COM = 1 class Leaf: def __init__(self): @@ -185,6 +199,7 @@ class Octree: def __init__(self): self.clear() self.min_aspect_ratio = 2.0/3.0 + self.type = Type.COV def clear(self): self.cell_boxes = [] @@ -226,14 +241,18 @@ def _progress(self, leaf_id, positions): self.leafs[leaf_id].id_first_child = self.cell_offsets[-1] - new_boxes = _create_combined_boxes(self.cell_boxes[leaf_id], self.min_aspect_ratio) + if self.type == Type.COV: + half = 0.5 * np.array([self.cell_boxes[leaf_id][i][1] + self.cell_boxes[leaf_id][i][0] for i in range(3)]) + elif self.type == Type.COM: + half = _get_centre_of_mass(self.permutations, positions, self.leafs[leaf_id].elem_offset, self.leafs[leaf_id].number_elements) + + new_boxes = _create_combined_boxes(self.cell_boxes[leaf_id], self.min_aspect_ratio, half) self.cell_offsets[-1] += len(new_boxes) self.leafs[leaf_id].number_children = len(new_boxes) for box in new_boxes: self.cell_boxes.append(box) - - #raise Exception() + else: pass diff --git a/setup.py b/setup.py index 800507d..60e7633 100644 --- a/setup.py +++ b/setup.py @@ -1,7 +1,7 @@ from setuptools import setup setup( name='dsmc', - version='0.9.0', + version='0.10.0', author='Leo Basov', python_requires='>=3.10, <4', packages=["dsmc", "dsmc.writer", "dsmc.mesh"], diff --git a/tests/unit/test_dsmc/octree.py b/tests/unit/test_dsmc/octree.py index 0fb249d..41d4843 100644 --- a/tests/unit/test_dsmc/octree.py +++ b/tests/unit/test_dsmc/octree.py @@ -223,29 +223,32 @@ def test__create_boxes(self): def test__get_min_aspect_ratio_1(self): box = np.array([(0.0, 1.0), (0.0, 10.0), (0.0, 100.0)]) + half = 0.5 * np.array([box[i][1] + box[i][0] for i in range(3)]) axis1 = 0 axis2 = 1 axis3 = 2 - self.assertEqual(0.01, oc._get_min_aspect_ratio(box, axis1)) - self.assertEqual(0.1, oc._get_min_aspect_ratio(box, axis2)) - self.assertEqual(10.0, oc._get_min_aspect_ratio(box, axis3)) + self.assertEqual(0.01, oc._get_min_aspect_ratio(box, axis1, half)) + self.assertEqual(0.1, oc._get_min_aspect_ratio(box, axis2, half)) + self.assertEqual(10.0, oc._get_min_aspect_ratio(box, axis3, half)) def test__get_min_aspect_ratio_2(self): box = np.array([(-1.0, 1.0), (-10.0, 10.0), (-100.0, 100.0)]) + half = 0.5 * np.array([box[i][1] + box[i][0] for i in range(3)]) axis1 = 0 axis2 = 1 axis3 = 2 - self.assertEqual(0.01, oc._get_min_aspect_ratio(box, axis1)) - self.assertEqual(0.1, oc._get_min_aspect_ratio(box, axis2)) - self.assertEqual(10.0, oc._get_min_aspect_ratio(box, axis3)) + self.assertEqual(0.01, oc._get_min_aspect_ratio(box, axis1, half)) + self.assertEqual(0.1, oc._get_min_aspect_ratio(box, axis2, half)) + self.assertEqual(10.0, oc._get_min_aspect_ratio(box, axis3, half)) def test__devide(self): box = np.array([(0.0, 1.0), (0.0, 10.0), (0.0, 100.0)]) - box_x1, box_x2 = oc._devide(box, 0) - box_y1, box_y2 = oc._devide(box, 1) - box_z1, box_z2 = oc._devide(box, 2) + half = 0.5 * np.array([box[i][1] + box[i][0] for i in range(3)]) + box_x1, box_x2 = oc._devide(box, 0, half) + box_y1, box_y2 = oc._devide(box, 1, half) + box_z1, box_z2 = oc._devide(box, 2, half) boxes = ((box_x1, box_x2), (box_y1, box_y2), (box_z1, box_z2)) @@ -269,9 +272,10 @@ def test__devide(self): def test__create_combined_boxes_1(self): box = np.array([(0.0, 1.0), (0.0, 10.0), (0.0, 100.0)]) + half = 0.5 * np.array([box[i][1] + box[i][0] for i in range(3)]) min_aspect_ratio = 0.0 boxes_old = oc._create_boxes(box) - boxes_new = oc._create_combined_boxes(box, min_aspect_ratio) + boxes_new = oc._create_combined_boxes(box, min_aspect_ratio, half) V = 0.0 for b in boxes_new: @@ -282,8 +286,9 @@ def test__create_combined_boxes_1(self): def test__create_combined_boxes_2(self): box = np.array([(0.0, 1.0), (0.0, 10.0), (0.0, 100.0)]) + half = 0.5 * np.array([box[i][1] + box[i][0] for i in range(3)]) min_aspect_ratio = 0.05 - boxes_new = oc._create_combined_boxes(box, min_aspect_ratio) + boxes_new = oc._create_combined_boxes(box, min_aspect_ratio, half) V = 0.0 for b in boxes_new: @@ -294,8 +299,9 @@ def test__create_combined_boxes_2(self): def test__create_combined_boxes_3(self): box = np.array([(-1.0, 1.0), (-10.0, 10.0), (-100.0, 100.0)]) + half = 0.5 * np.array([box[i][1] + box[i][0] for i in range(3)]) min_aspect_ratio = 0.05 - boxes_new = oc._create_combined_boxes(box, min_aspect_ratio) + boxes_new = oc._create_combined_boxes(box, min_aspect_ratio, half) V = 0.0 for b in boxes_new: @@ -306,8 +312,9 @@ def test__create_combined_boxes_3(self): def test__create_combined_boxes_4(self): box = np.array([(-1.0, 1.0), (-10.0, 10.0), (-100.0, 100.0)]) + half = 0.5 * np.array([box[i][1] + box[i][0] for i in range(3)]) min_aspect_ratio = 0.0 - boxes_new = oc._create_combined_boxes(box, min_aspect_ratio) + boxes_new = oc._create_combined_boxes(box, min_aspect_ratio, half) V = 0.0 for b in boxes_new: @@ -315,6 +322,19 @@ def test__create_combined_boxes_4(self): self.assertEqual(com.get_V(box), V) self.assertEqual(8, len(boxes_new)) + + + def test__get_centre_of_mass(self): + permutations = np.array([i for i in range(4)]) + positions = np.array([[-1.0, -1.0, 0.0], [1.0, -1.0, 0.0], [1.0, 1.0, 0.0], [-1.0, 1.0, 0.0]]) + offset = 0 + n_elements = 4 + + centre_of_mass = oc._get_centre_of_mass(permutations, positions, offset, n_elements) + + self.assertEqual(0.0, centre_of_mass[0]) + self.assertEqual(0.0, centre_of_mass[1]) + self.assertEqual(0.0, centre_of_mass[2]) class TestOctreeOctree(unittest.TestCase): def test_build(self):