|
| 1 | + |
| 2 | +__all__ = [ |
| 3 | + 'random_gaussian2d_model', |
| 4 | + 'random_unioned_circles_model', |
| 5 | + 'random_unioned_triangles_model', |
| 6 | +] |
| 7 | + |
| 8 | +from typing import Literal |
| 9 | +from collections import namedtuple |
| 10 | + |
| 11 | +from fealpy.backend import bm, TensorLike as Tensor |
| 12 | + |
| 13 | +from . import shapes as _S |
| 14 | + |
| 15 | + |
| 16 | +FunctionModel = namedtuple("FunctionModel", ["coef", "levelset", "shape"]) |
| 17 | + |
| 18 | + |
| 19 | +def _gaussian(p: Tensor, centers: Tensor, inv_cov: Tensor, base: float = 10.): |
| 20 | + # p (NC, NQ, GD) |
| 21 | + # centers (NGuassian, GD) |
| 22 | + # inv_cov (NGuassian, GD, GD) |
| 23 | + struct = p.shape[:-1] |
| 24 | + p = p.reshape(-1, p.shape[-1]) |
| 25 | + p0 = p[:, None, :] - centers[None, :, :] # (NC*NQ, NGuassian, GD) |
| 26 | + ind = bm.exp(-0.5 * bm.einsum("ngi, gij, ngj -> ng", p0, inv_cov, p0)) |
| 27 | + ind = bm.sum(ind, axis=-1).reshape(struct) # (NC, NQ) |
| 28 | + minim = bm.min(ind) |
| 29 | + maxim = bm.max(ind) |
| 30 | + ind = (ind - minim) / (maxim - minim) |
| 31 | + return base**ind |
| 32 | + |
| 33 | + |
| 34 | +def random_gaussian2d_model( |
| 35 | + num: int, |
| 36 | + box: tuple[float, float, float, float], |
| 37 | + major_lim: tuple[float, float], |
| 38 | + ecc_lim: tuple[float, float] = (0.5, 1.0), |
| 39 | + base: float = 10. |
| 40 | +): |
| 41 | + """ |
| 42 | + Generate a random Gaussian conductivity model. |
| 43 | +
|
| 44 | + Parameters: |
| 45 | + num (int): Number of Gaussian |
| 46 | + box (tuple[xmin, xmax, ymin, ymax]): The box of Gaussian centers |
| 47 | + major_lim (tuple[major_min, major_max]): The major axis covariance range |
| 48 | + ecc_lim (tuple[ecc_min, ecc_max]): The eccentricity range |
| 49 | + base (float): The ratio of the conductivity value that ranges [1.0, ratio] |
| 50 | +
|
| 51 | + Returns: |
| 52 | + namedtuple |
| 53 | + - coef (Callable): A function that takes points and returns the conductivity values. |
| 54 | + - levelset (None): No levelset function. |
| 55 | + """ |
| 56 | + gaussian = _S.random_gaussian_2d( |
| 57 | + num, box, major_lim, ecc_lim |
| 58 | + ) |
| 59 | + |
| 60 | + def gaussian_conductivity(points: Tensor, *args, **kwargs) -> Tensor: |
| 61 | + return _gaussian( |
| 62 | + points, |
| 63 | + bm.from_numpy(gaussian.centers), |
| 64 | + bm.from_numpy(gaussian.invcov), |
| 65 | + base |
| 66 | + ) |
| 67 | + |
| 68 | + gaussian_conductivity.__dict__['coordtype'] = 'cartesian' |
| 69 | + |
| 70 | + return FunctionModel(gaussian_conductivity, None, gaussian) |
| 71 | + |
| 72 | + |
| 73 | +def random_unioned_circles_model( |
| 74 | + num: int, |
| 75 | + values: tuple[float, float] = (10.0, 1.0), |
| 76 | +): |
| 77 | + """Generate a random unioned circles conductivity model. |
| 78 | +
|
| 79 | + Parameters: |
| 80 | + num (int): Number of circles |
| 81 | + values (tuple[float, float]): The conductivity values of the circles and the rest. |
| 82 | +
|
| 83 | + Returns: |
| 84 | + namedtuple: |
| 85 | + - coef (Callable): A function that takes points and returns the conductivity values. |
| 86 | + - levelset (Callable): Level set function. |
| 87 | + """ |
| 88 | + circles = _S.random_circles(num) |
| 89 | + level_set_func = lambda p: _S.circle_union_levelset(p, circles) |
| 90 | + |
| 91 | + def circle_conductivity(points: Tensor, *args, **kwargs) -> Tensor: |
| 92 | + inclusion = level_set_func(points) < 0. # a bool tensor on quadrature points. |
| 93 | + sigma = bm.empty(points.shape[:2], **bm.context(points)) # (Q, C) |
| 94 | + sigma = bm.set_at(sigma, inclusion, values[0]) |
| 95 | + sigma = bm.set_at(sigma, ~inclusion, values[1]) |
| 96 | + return sigma |
| 97 | + |
| 98 | + circle_conductivity.__dict__['coordtype'] = 'cartesian' |
| 99 | + |
| 100 | + return FunctionModel(circle_conductivity, level_set_func, circles) |
| 101 | + |
| 102 | + |
| 103 | +def random_unioned_triangles_model( |
| 104 | + num: int, |
| 105 | + box: tuple[float, float, float, float] = [-1., 1., -1., 1.], |
| 106 | + kind: Literal["", "equ"] = "", |
| 107 | + rlim: tuple[float, float] = (0.5, 1.0), |
| 108 | + values: tuple[float, float] = (10.0, 1.0), |
| 109 | +): |
| 110 | + """Generate a random unioned triangles conductivity model. |
| 111 | +
|
| 112 | + Parameters: |
| 113 | + num (int): Number of triangles |
| 114 | + box (tuple[xmin, xmax, ymin, ymax]): The box of triangles |
| 115 | + kind (str): The kind of triangles. |
| 116 | + rlim (tuple[r_min, r_max]): The radius range, available when kind is "equ". |
| 117 | + values (tuple[float, float]): The conductivity values of the triangles and the rest. |
| 118 | +
|
| 119 | + Returns: |
| 120 | + namedtuple: |
| 121 | + - coef (Callable): A function that takes points and returns the conductivity values. |
| 122 | + - levelset (Callable): Level set function. |
| 123 | + """ |
| 124 | + if kind == "equ": |
| 125 | + tris = _S.random_equ_triangles(num, box[0:2], box[2:4], rlim) |
| 126 | + else: |
| 127 | + tris = _S.random_ccw_triangles(num, box[0:2], box[2:4]) |
| 128 | + level_set_func = lambda p: _S.triangle_union_levelset(p, tris) |
| 129 | + |
| 130 | + def triangle_conductivity(points: Tensor, *args, **kwargs) -> Tensor: |
| 131 | + inclusion = level_set_func(points) < 0. # a bool tensor on quadrature points. |
| 132 | + sigma = bm.empty(points.shape[:2], **bm.context(points)) # (Q, C) |
| 133 | + sigma = bm.set_at(sigma, inclusion, values[0]) |
| 134 | + sigma = bm.set_at(sigma, ~inclusion, values[1]) |
| 135 | + return sigma |
| 136 | + |
| 137 | + triangle_conductivity.__dict__['coordtype'] = 'cartesian' |
| 138 | + |
| 139 | + return FunctionModel(triangle_conductivity, level_set_func, tris) |
0 commit comments