From 9b7457196f044aa7de35f83fb78fc339d1b615b4 Mon Sep 17 00:00:00 2001 From: Jakob Date: Tue, 9 Jan 2024 16:11:15 +0100 Subject: [PATCH] rm dir; start over --- src/data/helmholtz/meshing/Domain.py | 158 ------------------ src/data/helmholtz/meshing/DomainManager.py | 10 -- src/data/helmholtz/meshing/MeshStrategy.py | 30 ---- src/data/helmholtz/meshing/Property.py | 155 ----------------- src/data/helmholtz/meshing/shapes/CShape.py | 54 ------ .../meshing/shapes/CShapeUnitCell.py | 66 -------- .../helmholtz/meshing/shapes/Rectangle.py | 64 ------- src/data/helmholtz/meshing/shapes/Shape.py | 82 --------- src/data/helmholtz/meshing/shapes/__init__.py | 0 src/data/helmholtz/meshing/util.py | 15 -- 10 files changed, 634 deletions(-) delete mode 100644 src/data/helmholtz/meshing/Domain.py delete mode 100644 src/data/helmholtz/meshing/DomainManager.py delete mode 100644 src/data/helmholtz/meshing/MeshStrategy.py delete mode 100644 src/data/helmholtz/meshing/Property.py delete mode 100644 src/data/helmholtz/meshing/shapes/CShape.py delete mode 100644 src/data/helmholtz/meshing/shapes/CShapeUnitCell.py delete mode 100644 src/data/helmholtz/meshing/shapes/Rectangle.py delete mode 100644 src/data/helmholtz/meshing/shapes/Shape.py delete mode 100644 src/data/helmholtz/meshing/shapes/__init__.py delete mode 100644 src/data/helmholtz/meshing/util.py diff --git a/src/data/helmholtz/meshing/Domain.py b/src/data/helmholtz/meshing/Domain.py deleted file mode 100644 index dc5072db..00000000 --- a/src/data/helmholtz/meshing/Domain.py +++ /dev/null @@ -1,158 +0,0 @@ -from __future__ import annotations - -from collections import defaultdict -from typing import List, Optional, Tuple - -import numpy as np - -from .Property import AdiabaticAbsorberProperty, Property -from .util import Direction - - -class BoxDomain: - """Maintains a subspace of the domain (which may be the entire domain) - - Top most domain will be placed at (0, 0) with its lower left corner. This class also manages its subdomains. - The locations of the subdomains are placed in a relative context to this one. - - """ - - def __init__( - self, - name: str = "BoxDomain", - dx: Optional[float] = None, - dy: Optional[float] = None, - properties: dict[str, Property] = None, - sub_domains: Optional[List[BoxDomain]] = None, - direction: Optional[Direction] = Direction.positive_x, - ): - """ - - Args: - name: to identify domain - dx: extend in positive x direction - dy: extend in positive y direction - properties: dict, name of property with a property object, Properties defined on this domain - sub_domains: List, of domains encompassed by this domain - direction: Direction, domains are stacked as they are listed in sub_domains in this direction - """ - if not any([properties, all([dx, dy])]): - # if neither is provided, this would be an empty object - raise ValueError("Domain is not well-defined.") - self.name = name - self.dimension = [dx, dy] - if not all(self.dimension): - self.dimension = [0.0, 0.0] - if properties is None: - properties = defaultdict(lambda: None) - self.properties = properties # holds information of physical parameters within the domain - if sub_domains: - self.sub_domains = [ - [ - s, - [ - 0.0, - ] - * 2, - ] - for s in sub_domains - ] # (absolute bbox coordinates in domain, sub_domain) - else: - self.sub_domains = [] - self.direction = (direction + 1) % 2 - self.orthogonal_direction = direction % 2 - - # set properties - self.update_coordinates() - - # bounding box - self.bbox = [0.0, 0.0] # coordinate of upper right corner - if self.sub_domains: - self.bbox = self.sub_domains[-1][0].dimension.copy() # orthogonal dimension - self.bbox[self.direction] += self.sub_domains[-1][1][self.direction] # sum of dx and position - else: - self.bbox = self.dimension - - def update_coordinates(self) -> None: - """Sets coordinates within the contex of all domains. - - Domains are stacked in the following order: 1) own domain 2) sub-domains in direction in the order of the list - in self.sub_domains. - - Returns: - - """ - coordinates = [ - 0.0, - ] * 2 - - # the dimension of this object is the first. - coordinates[self.direction] += self.dimension[self.direction] - - # shift children accordingly - if not self.sub_domains: - return - # check data alignment - dimensions = set([self.dimension[self.orthogonal_direction]]) - for domain, _ in self.sub_domains: - dimensions.add(domain.dimension[self.orthogonal_direction]) - dimensions.discard(0.0) # domain without own "body" - if len(dimensions) > 1: - raise ValueError("Only rectangular domains allowed!") - - for i, (domain, _) in enumerate(self.sub_domains): - # modify sub_domain information to match relative position from this domain - self.sub_domains[i][1] = coordinates.copy() # points at the lower left corner - coordinates[self.direction] += domain.bbox[self.direction] - - def eval(self, property_name: str, x) -> np.array: - """ - - Args: - property_name: physical property definition (as a domain can hold multiple different ones) - x: tuple, position - - Returns: - - """ - output_value = np.zeros(x[0].shape, dtype=np.complex128) - for sub_domain, pos in self.sub_domains: - relative_x = [x[0] - pos[0], x[1] - pos[1]] - # only values inside - output_value += sub_domain.eval(property_name, relative_x) - if self.properties[property_name]: - inside = (0.0 <= x[0]) * (x[0] <= self.bbox[0]) * (0.0 <= x[1]) * (x[1] <= self.bbox[1]) - prop = self.properties[property_name] - output_value += prop.eval(x) * inside - return output_value - - -class AdiabaticLayer(BoxDomain): - """Layer of an Adiabatic Absorber - - Adiabatic Absrobers are an extension to perfectly matched layers. They were introduced as PMLs sometimes fail to be - reflectionless even in the limit of infinite resolution. These absorbers make the reflections negligible by - gradually increasing the material absorption. - Details can be found in Oskooi, A. F., Zhang, L., Avniel, Y. & Johnson, S. G. The failure of perfectly matched - layers, and towards their redemption by adiabatic absorbers. Opt. Express 16, 11376–11392 (2008). - """ - - def __init__( - self, - dx: float, - dy: float, - round_trip: float, - direction_depth: dict[int, List[float, float]], - properties: List[Tuple[str, float]], - ): - """ - - Args: - dx: extend in x direction - dy: extend in y direction - round_trip: float, governs the size of reflections caused by this layer - direction_depth: dict, key: direction identifier, value: Starting and end point (key needs to be positive) - properties: List of Tuples, defines which quantities/Properties should be modified by this domain - """ - props = {prop[0]: AdiabaticAbsorberProperty(2, round_trip, prop[1], direction_depth) for prop in properties} - super().__init__("Adiabaitc Absorber", dx, dy, props) diff --git a/src/data/helmholtz/meshing/DomainManager.py b/src/data/helmholtz/meshing/DomainManager.py deleted file mode 100644 index a20eb867..00000000 --- a/src/data/helmholtz/meshing/DomainManager.py +++ /dev/null @@ -1,10 +0,0 @@ -class MeshManager: - """manages domain descriptions and meshes""" - - def __init__(self): - self.mesh = None - self.boundaries = None - self.mesh_strategy = None - - def get_mesh(self): - pass diff --git a/src/data/helmholtz/meshing/MeshStrategy.py b/src/data/helmholtz/meshing/MeshStrategy.py deleted file mode 100644 index f16e394c..00000000 --- a/src/data/helmholtz/meshing/MeshStrategy.py +++ /dev/null @@ -1,30 +0,0 @@ -import pathlib -from abc import ABC, abstractmethod - -import dolfinx - -from .Domain import BoxDomain - - -class MeshStrategy(ABC): - @abstractmethod - def load_mesh(self) -> dolfinx.mesh.Mesh: - pass - - -class MeshFromFile(MeshStrategy): - def __init__(self, data_file: pathlib.Path): - super().__init__() - - def load_mesh(self) -> dolfinx.mesh.Mesh: - pass - - -class MeshFromBox(MeshStrategy): - """creates a mesh within a box like dolfinx.mesh.Mesh structure""" - - def __init__(self, domain: BoxDomain): - self.domain = domain - - def load_mesh(self) -> dolfinx.mesh.Mesh: - pass diff --git a/src/data/helmholtz/meshing/Property.py b/src/data/helmholtz/meshing/Property.py deleted file mode 100644 index 59a394cf..00000000 --- a/src/data/helmholtz/meshing/Property.py +++ /dev/null @@ -1,155 +0,0 @@ -from abc import ABC, abstractmethod -from collections import defaultdict -from typing import List - -import numpy as np - - -class Property(ABC): - """Defines a property on a domain - - Properties modify physical quantities of domains. These may be structures of physical parameters. - - """ - - @abstractmethod - def eval(self, x) -> np.array: - """ - - Args: - x: tuple, position - - Returns: value by which property should be modified at location x - - """ - pass - - -class ConstantProperty(Property): - """Constant. - - Returns always the same constant value. - - """ - - def __init__(self, const: float): - """ - - Args: - const: value of this property - """ - self.value = const - - def eval(self, x) -> np.array: - """ - - Args: - x: location vector - - Returns: the value of this property at the position x - - """ - return self.value * np.ones(x[0].shape) - - -class CylindricalCrystalProperty(Property): - """Equally spaced crystals within a given rectangle - - Equally spaced cylindrical crystals within a rectangular domain. - They are equally spaced to each other and to the boundaries of the rectangle. - The cylinders manifest in a modification (by addition) of a property defined on the domain. - Overlapping crystals are treated as a union of both (i.e., they do not amplify that region further). - """ - - def __init__( - self, - dx: float, - dy: float, - n_x: int, - n_y: int, - radius: float, - value: float, - ): - """ - - Args: - dx: size of rectangle in x direction - dy: size of rectangle in y direction - n_x: number of crystals in x direction - n_y: number of crystals in y direction - radius: radius of each crystal - value: modification (by addition) of a given property within domain - """ - stride = (dx / n_x, dy / n_y) # start in rect (if it fits) - self.centers = [] - for row in range(n_x): - for col in range(n_y): - self.centers.append((stride[0] * (row + 0.5), stride[1] * (col + 0.5))) - self.squared_radius = radius**2 - self.value = value - - def eval(self, x) -> np.array: - """ - - Args: - x: tuple/vector, position - - Returns: Value of the modified property at location x - - """ - in_circle = np.max( - np.array([(x[0] - c[0]) ** 2 + (x[1] - c[1]) ** 2 <= self.squared_radius for c in self.centers]), - axis=0, - ) - return in_circle * self.value - - -class AdiabaticAbsorberProperty(Property): - """Modifies wavenumber in property domain to truncate wave - - Adiabatic Absrobers are an extension to perfectly matched layers. They were introduced as PMLs sometimes fail to be - reflectionless even in the limit of infinite resolution. These absorbers make the reflections negligible by - gradually increasing the material absorption. - Details can be found in Oskooi, A. F., Zhang, L., Avniel, Y. & Johnson, S. G. The failure of perfectly matched - layers, and towards their redemption by adiabatic absorbers. Opt. Express 16, 11376–11392 (2008). - """ - - def __init__( - self, - degree: int, - round_trip: float, - value: float, - direction_depth: dict[int, List[float]], - ): - self.direction_properties = defaultdict(lambda: [-1.0, -1.0]) - self.direction_properties.update(direction_depth) # key: direction, value: [start, end] - self.degree = degree - self.round_trip = round_trip - self.sigma_0 = { - axis: -(self.degree + 1) * np.log(self.round_trip) / (2.0 * (pos[1] - pos[0])) * (-1) ** (pos[1] > pos[0]) - for axis, pos in self.direction_properties.items() - } - self.value = value - - def eval(self, x) -> np.array: - """ - - Args: - x: tuple/vector, location on which the property is evaluated - - Returns: Modified property value on location x - - """ - xi = { - direction: (x[direction.get_axis()] - pos[0]) / (pos[1] - pos[0]) - for direction, pos in self.direction_properties.items() - } - # bound xi as values above 1 and below 0 indicate that xi is not in absorption range - for direction, value in xi.items(): - inside = (value >= 0) * (value <= 1) - xi[direction] = value * inside - sigma_x = [] - for direction, zeta in xi.items(): - sigma_x.append(self.sigma_0[direction] * (zeta**self.degree)) - sigma_x = np.sum(np.array(sigma_x), axis=0) - return 2j * sigma_x * self.value - sigma_x**2 diff --git a/src/data/helmholtz/meshing/shapes/CShape.py b/src/data/helmholtz/meshing/shapes/CShape.py deleted file mode 100644 index 2113e823..00000000 --- a/src/data/helmholtz/meshing/shapes/CShape.py +++ /dev/null @@ -1,54 +0,0 @@ -from operator import add - -from src.utility import run_once - -from .Shape import Shape - - -class CShape(Shape): - """2D C-shape - - C-shape is a specific kind of sonic crystal. They are usually arranged in a periodic fashion. A grid of these - c-shapes can exhibit sub-frequency behavior, specific to meta-materials. An overview over meta-materials can be - found in "Acoustic meta-materials: From local resonances to broad horizons" by Ma et al. - - Attributes: - - """ - - def __init__(self, x: float, y: float, z: float, r_outer: float, r_inner: float, b: float): - """ - - Args: - x: x-coordinate of center (defined as center of both arcs) - y: y-coordinate of center (defined as center of both arcs) - z: z-coordinate of center (defined as center of both arcs) - r_outer: outer radius - r_inner: inner radius - b: width of the slit - """ - super().__init__(x, y, z, r_outer, r_inner, b) - - @run_once - def generate(self, x: float, y: float, z: float, r_outer: float, r_inner: float, b: float): - """ - Args: - x: x-coordinate of center (defined as center of both arcs) - y: y-coordinate of center (defined as center of both arcs) - z: z-coordinate of center (defined as center of both arcs) - r_outer: outer radius - r_inner: inner radius - b: width of the slit - """ - coordinates = [x, y, z] - r = sorted([r_inner, r_outer]) - rect_offset = [0.0, -b / 2.0, 0.0] - - # basic shapes - outer_disk = self.factory.addDisk(*coordinates, r[1], r[1]) - inner_disk = self.factory.addDisk(*coordinates, r[0], r[0]) - rectangle_coordinates = list(map(add, coordinates, rect_offset)) - slit = self.factory.addRectangle(*rectangle_coordinates, -r[1], b) - - # cut inner circle to get pipe and slit pipe to obtain C-shape - self.factory.cut([(2, outer_disk)], [(2, inner_disk), (2, slit)]) diff --git a/src/data/helmholtz/meshing/shapes/CShapeUnitCell.py b/src/data/helmholtz/meshing/shapes/CShapeUnitCell.py deleted file mode 100644 index e3785658..00000000 --- a/src/data/helmholtz/meshing/shapes/CShapeUnitCell.py +++ /dev/null @@ -1,66 +0,0 @@ -from src.utility import run_once - -from .CShape import CShape -from .Rectangle import Rectangle -from .Shape import Shape - - -class CShapeUnitCell(Shape): - """Unit Cell or crystal with a c-shape inside - - Attributes: - - """ - - def __init__( - self, x: float, y: float, z: float, dx: float, dy: float, r_outer: float, r_inner: float, b, cut: bool = False - ): - """Unit cell with c-shape placed at the center - - Args: - x: coordinate - y: coordinate - z: coordinate - dx: x-dimension - dy: y-dimension - r_outer: c-shape outer radius - r_inner: c-shape inner radius - b: slit width of c-shape - cut: c-shape should be cut from cell mesh - """ - raise DeprecationWarning("DO not use") - super().__init__(x, y, z, dx, dy, r_outer, r_inner, b, cut) - - @run_once - def generate( - self, x: float, y: float, z: float, dx: float, dy: float, r_outer: float, r_inner: float, b, cut: bool = True - ) -> None: - """ - - Args: - x: coordinate - y: coordinate - z: coordinate - dx: x-dimension - dy: y-dimension - r_outer: c-shape outer radius - r_inner: c-shape inner radius - b: slit width of c-shape - cut: c-shape should be cut from cell mesh - - Returns: None - - """ - center_x = x + dx / 2 - center_y = y + dy / 2 - - domain = Rectangle(x, y, z, dx, dy) - c_shape = CShape(center_x, center_y, z, r_outer, r_inner, b) - - # cut or fragment for final domain - self.children[domain.__class__.__name__] = domain - if cut: - self.factory.cut(domain(), c_shape()) - return - self.factory.fragment(domain(), c_shape()) - self.children[c_shape.__class__.__name__] = c_shape diff --git a/src/data/helmholtz/meshing/shapes/Rectangle.py b/src/data/helmholtz/meshing/shapes/Rectangle.py deleted file mode 100644 index 5f076627..00000000 --- a/src/data/helmholtz/meshing/shapes/Rectangle.py +++ /dev/null @@ -1,64 +0,0 @@ -import gmsh -import numpy as np - -from src.utility import run_once - -from .Shape import Shape - - -class Rectangle(Shape): - """2D Rectangle - - Attributes: - - """ - - def __init__(self, x: float, y: float, z: float, dx: float, dy: float): - """ - - Args: - x: lower left corner coordinate - y: lower left corner coordinate - z: lower left corner coordinate - dx: x-size - dy: y-size - """ - super().__init__(x, y, z, dx, dy) - - # boundary tags ar not always assigned in the same order - self.boundaries = { - "top": None, - "left": None, - "bottom": None, - "right": None, - } - a_tol = min(dx, dy) / 1e5 - for boundary in self.elements[1]: - coordinates = gmsh.model.getBoundingBox(1, boundary) - delta_x = coordinates[3] - coordinates[0] - - if np.isclose(delta_x, dx, atol=a_tol): - # top or bottom - if np.isclose(coordinates[1], y, atol=a_tol): - self.boundaries["bottom"] = boundary - continue - self.boundaries["top"] = boundary - continue - # left or right - if np.isclose(coordinates[0], x, atol=a_tol): - self.boundaries["left"] = boundary - continue - self.boundaries["right"] = boundary - - @run_once - def generate(self, x: float, y: float, z: float, dx: float, dy: float): - """ - - Args: - x: lower left corner coordinate - y: lower left corner coordinate - z: lower left corner coordinate - dx: x-size - dy: y-size - """ - self.factory.addRectangle(x, y, z, dx, dy) diff --git a/src/data/helmholtz/meshing/shapes/Shape.py b/src/data/helmholtz/meshing/shapes/Shape.py deleted file mode 100644 index 6e2d3fcd..00000000 --- a/src/data/helmholtz/meshing/shapes/Shape.py +++ /dev/null @@ -1,82 +0,0 @@ -import warnings -from abc import ABC, abstractmethod -from collections import defaultdict - -import gmsh - -from src.utility import run_once - - -class Shape(ABC): - """Shape for gmsh meshing and cutting - - Attributes: - - """ - - def __init__(self, *args, **kwargs): - self.factory = gmsh.model.occ # it is not encouraged to mix interfaces - self.elements = defaultdict(lambda: set()) # holds elements associated with the top level shape (keys are dims) - self.children = {} - - # generate shape - gmsh.model.occ.synchronize() - elements_before = gmsh.model.get_entities(-1) - self.generate(*args, **kwargs) - gmsh.model.occ.synchronize() - elements_after = gmsh.model.get_entities(-1) - new_elements = set(elements_after) - set(elements_before) - - # update elements - for element in new_elements: - self.elements[element[0]].add(element[1]) - - # child elements need to be accessed via child - for dim in self.elements: - for child in self.children.values(): - self.elements[dim] -= child.elements[dim] - - @abstractmethod - @run_once - def generate(self, *args, **kwargs) -> None: - """Initializes the entire shape, build all children and populates attributes of shape - - Args: - *args: - **kwargs: - - Returns: - - """ - - def __call__(self, dim: int = None): - """Get dimensions and tags identifying this shape - - Returns: List of tuples (dim, tag) pertaining to this shape (children excluded) - - """ - if not self.elements: - # this shape holds no elements itself - collect recursively from children - if not self.children: - warnings.warn(f"Shape {self.__class__.__name__} holds not elements, but was called.") - return [] - out = [] - for child in self.children.values(): - out += child() - return out - if not dim: - # the highest dimensional element usually defines a shape - dim = max(self.elements) - return [(dim, element) for element in self.elements[dim]] - - def add_physical_group(self, dim: int, name: str) -> None: - """Add physical group name to all elements of dim - - Args: - dim: for elements of this dim - name: physical group - - Returns: - - """ - gmsh.model.add_physical_group(dim, list(self.elements[dim]), name=name) diff --git a/src/data/helmholtz/meshing/shapes/__init__.py b/src/data/helmholtz/meshing/shapes/__init__.py deleted file mode 100644 index e69de29b..00000000 diff --git a/src/data/helmholtz/meshing/util.py b/src/data/helmholtz/meshing/util.py deleted file mode 100644 index 20c1fcee..00000000 --- a/src/data/helmholtz/meshing/util.py +++ /dev/null @@ -1,15 +0,0 @@ -from enum import IntEnum - - -class Direction(IntEnum): - """Direction in a 3-dimensional space.""" - - positive_x = 1 # start at 1 to utilize symmetry in both directions - positive_y = 2 - positive_z = 3 - negative_x = -1 - negative_y = -2 - negative_z = -3 - - def get_axis(self): - return abs(self.value) - 1