From 892231d3aac56006303b105b48ceea8aae922e20 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Thu, 24 Aug 2023 15:52:37 +0200 Subject: [PATCH 1/5] add cs module, mesh.csgmsh function --- nutils/cs/__init__.py | 65 ++++++ nutils/cs/_axes.py | 134 +++++++++++ nutils/cs/_field.py | 158 +++++++++++++ nutils/cs/_gmsh_stub.py | 53 +++++ nutils/cs/_shape.py | 485 ++++++++++++++++++++++++++++++++++++++++ nutils/cs/_util.py | 15 ++ nutils/mesh.py | 18 +- tests/test_cs.py | 406 +++++++++++++++++++++++++++++++++ 8 files changed, 1333 insertions(+), 1 deletion(-) create mode 100644 nutils/cs/__init__.py create mode 100644 nutils/cs/_axes.py create mode 100644 nutils/cs/_field.py create mode 100644 nutils/cs/_gmsh_stub.py create mode 100644 nutils/cs/_shape.py create mode 100644 nutils/cs/_util.py create mode 100644 tests/test_cs.py diff --git a/nutils/cs/__init__.py b/nutils/cs/__init__.py new file mode 100644 index 000000000..f4b156967 --- /dev/null +++ b/nutils/cs/__init__.py @@ -0,0 +1,65 @@ +'''Constructive Solid Geometry helper for Gmsh''' + + +from ._shape import ( + Box, + Circle, + Cylinder, + Ellipse, + Entities, + Interval, + Line, + Path, + Point, + Rectangle, + Skeleton, + Sphere, + generate_mesh, +) + + +from ._field import ( + AsField, + Ball, + LocalRefinement, + Max, + Min, + set_background, + x, + y, + z, +) + + +def write(fname: str, entities: Entities, elemsize: AsField, order: int = 1) -> None: + 'Create .msh file based on Constructive Solid Geometry description.' + + import treelog, gmsh # type: ignore + + gmsh.initialize() + try: + gmsh.option.setNumber('General.Terminal', 0) + gmsh.option.setNumber('Mesh.Binary', 1) + gmsh.option.setNumber('Mesh.ElementOrder', order) + gmsh.option.setNumber('Mesh.CharacteristicLengthExtendFromBoundary', 0) + gmsh.option.setNumber('Mesh.CharacteristicLengthFromPoints', 0) + gmsh.option.setNumber('Mesh.CharacteristicLengthFromCurvature', 0) + + gmsh.logger.start() + try: + set_background(gmsh.model.mesh.field, elemsize) + generate_mesh(gmsh.model, entities) + gmsh.write(fname) + finally: + with treelog.context('gmsh'): + for line in gmsh.logger.get(): + level, sep, msg = line.partition(': ') + if level in ('Debug', 'Info', 'Warning', 'Error'): + getattr(treelog, level.lower())(msg) + gmsh.logger.stop() + + finally: + gmsh.finalize() + + +# vim:sw=4:sts=4:et diff --git a/nutils/cs/_axes.py b/nutils/cs/_axes.py new file mode 100644 index 000000000..47e825fae --- /dev/null +++ b/nutils/cs/_axes.py @@ -0,0 +1,134 @@ +from ._gmsh_stub import OCC, DimTags +from typing import Tuple, Any, Union +import numpy +import math + +XY = Tuple[float,float] +XYZ = Tuple[float,float,float] + + +class Axes2: + 'representation of 2D, positive, orthonormal axes' + + @classmethod + def eye(cls) -> 'Axes2': + return cls(0.) + + @classmethod + def from_x(cls, xaxis: XY) -> 'Axes2': + cosθ, sinθ = xaxis + return cls(math.atan2(sinθ, cosθ)) + + def __init__(self, rotation: float): + self._θ = rotation + + def __len__(self) -> int: + return 2 + + def __getitem__(self, s): + sin = numpy.sin(self._θ) + cos = numpy.cos(self._θ) + return numpy.array([[cos, sin], [-sin, cos]][s]) + + def rotate(self, angle: float) -> 'Axes2': + return Axes2(self._θ + angle) + + def as_3(self, origin: XY) -> Tuple['Axes3', XYZ]: + return Axes3.from_rotation_vector((0, 0, self._θ)), (*origin, 0) + + +class Axes3: + 'representation of 3D, positive, orthonormal axes' + # immutable representation of a unit quaternion, see https://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation + + @classmethod + def eye(cls) -> 'Axes3': + return cls(1., (0., 0., 0.)) + + @classmethod + def from_rotation_vector(cls, v: XYZ) -> 'Axes3': + v_ = numpy.asarray(v) + θ = numpy.linalg.norm(v_) / 2 + return cls(numpy.cos(θ), v_ * (-.5 * numpy.sinc(θ/numpy.pi))) + + @classmethod + def from_xy(cls, xaxis: XYZ, yaxis: XYZ) -> 'Axes3': + Q = numpy.array([xaxis, yaxis, numpy.cross(xaxis, yaxis)]) + + K = numpy.empty([4,4]) + K[1:,1:] = Q + Q.T + K[0,1:] = K[1:,0] = Q[[2,0,1],[1,2,0]] - Q[[1,2,0],[2,0,1]] + K[0,0] = 2 * numpy.trace(Q) + + eigval, eigvec = [v.T[-1] for v in numpy.linalg.eigh(K)] + # is xaxis and yaxis are orthonormal then (eigval-K[0,0]/2)/3 == 1 + + return cls(eigvec[0], eigvec[1:]) + + def __init__(self, w: float, v: XYZ): + self._w = w + self._v = numpy.array(v) + + def __len__(self) -> int: + return 3 + + def __getitem__(self, s): + I = numpy.eye(3)[s] + return 2 * ((.5 - self._v @ self._v) * I + self._v[s, numpy.newaxis] * self._v + numpy.cross(I, self._w * self._v)) + + def rotate(self, rotvec: XYZ) -> 'Axes3': + other = Axes3.from_rotation_vector(rotvec) + # hamilton product + return Axes3(self._w * other._w - self._v @ other._v, self._w * other._v + other._w * self._v + numpy.cross(self._v, other._v)) + + def as_3(self, origin: XYZ) -> Tuple['Axes3', XYZ]: + return self, origin + + @property + def rotation_axis(self) -> XYZ: + return tuple(self._v) + + @property + def rotation_angle(self) -> float: + return -2 * math.atan2(numpy.linalg.norm(self._v), self._w) + + @property + def rotation_vector(self) -> XYZ: + n = numpy.linalg.norm(self._v) + return self._v * (n and -2 * math.atan2(n, self._w) / n) + + +class Orientation: + + def __init__(self, origin: Union[XY, XYZ], **kwargs: Any): + ndims = len(origin) + args = ', '.join(sorted(kwargs)) + if args == 'axes': + axes = kwargs['axes'] + assert len(axes) == ndims + elif ndims == 2 and not args: + axes = Axes2.eye() + elif ndims == 3 and not args: + axes = Axes3.eye() + elif ndims == 2 and args == 'rotation': + axes = Axes2(kwargs['rotation']) + elif ndims == 2 and args == 'xaxis': + axes = Axes2.from_x(kwargs['xaxis']) + elif ndims == 3 and args == 'xaxis, yaxis': + axes = Axes3.from_xy(kwargs['xaxis'], kwargs['yaxis']) + elif ndims == 3 and args == 'rotvec': + axes = Axes3.from_rotation_vector(kwargs['rotvec']) + else: + raise ValueError(f'cannot create {ndims}D orientation based on arguments {args}') + self.origin = numpy.array(origin) + self.ndims = ndims + self.axes = axes + + def orient(self, occ: OCC, dimtags: DimTags) -> None: + 'position a shape in 3D space' + + axes, origin = self.axes.as_3(self.origin) + if any(origin): + occ.translate(dimtags, *origin) + if axes.rotation_angle: + occ.rotate(dimtags, *origin, *axes.rotation_axis, -axes.rotation_angle) diff --git a/nutils/cs/_field.py b/nutils/cs/_field.py new file mode 100644 index 000000000..07392a3e9 --- /dev/null +++ b/nutils/cs/_field.py @@ -0,0 +1,158 @@ +from typing import Tuple, Union, Protocol, runtime_checkable, NewType +from dataclasses import dataclass + + +# STUB FOR RELEVANT GMSH API + +Tag = NewType('Tag', int) + +class FieldFactory(Protocol): + def add(self, name: str) -> Tag: + ... + def setString(self, tag: Tag, option: str, value: str) -> None: + ... + def setNumber(self, tag: Tag, option: str, value: Union[int, float]) -> None: + ... + def setAsBackgroundMesh(self, tag: Tag) -> None: + ... + +# END STUB + + +AsField = Union['Field',int,float] + + +@runtime_checkable +class Field(Protocol): + 'A scalar field over a coordinate system.' + + def getexpr(self, ff: FieldFactory) -> str: + ... + def __add__(self, other: AsField) -> 'Field': + return FieldOp(self, as_field(other), '+') + def __radd__(self, other: AsField) -> 'Field': + return FieldOp(as_field(other), self, '+') + def __sub__(self, other: AsField) -> 'Field': + return FieldOp(self, as_field(other), '-') + def __rsub__(self, other: AsField) -> 'Field': + return FieldOp(as_field(other), self, '-') + def __mul__(self, other: AsField) -> 'Field': + return FieldOp(self, as_field(other), '*') + def __rmul__(self, other: AsField) -> 'Field': + return FieldOp(as_field(other), self, '*') + def __truediv__(self, other: AsField) -> 'Field': + return FieldOp(self, as_field(other), '/') + + +def as_field(f: AsField) -> Field: + if isinstance(f, Field): + return f + if isinstance(f, (int, float)): + return Constant(f) + raise ValueError(f'cannot interpret {f!r} as a field') + + +@dataclass(frozen=True) +class FieldOp(Field): + + a: Field + b: Field + op: str + + def getexpr(self, ff: FieldFactory) -> str: + a = self.a.getexpr(ff) + b = self.b.getexpr(ff) + return f'({a}){self.op}({b})' + + +@dataclass(frozen=True) +class Constant(Field): + 'Constant element size' + + size: float + + def getexpr(self, ff: FieldFactory) -> str: + return str(self.size) + + +@dataclass(frozen=True) +class Coord(Field): + + coord: str + + def getexpr(self, ff: FieldFactory) -> str: + return self.coord + +x = Coord('x') +y = Coord('y') +z = Coord('z') + + +@dataclass(frozen=True) +class LocalRefinement(Field): + 'Refine elements according to a gaussian bell.' + + center: Tuple[float,float] = (0., 0.) + radius: float = 1. + factor: float = 2. + + def getexpr(self, ff: FieldFactory) -> str: + x, y = self.center + s = 1 - 1/self.factor + return f'1-{s}*exp(-((x-{x})^2+(y-{y})^2)/({self.radius})^2)' + + +@dataclass(frozen=True) +class Ball(Field): + 'Refine elements uniformly inside a circle' + + center: Tuple[float,float] + radius: float + inside: float + outside: float + thickness: float = 0. + + def getexpr(self, ff: FieldFactory) -> str: + tag = ff.add('Ball') + ff.setNumber(tag, 'Radius', self.radius) + ff.setNumber(tag, 'Thickness', self.thickness) + ff.setNumber(tag, 'VIn', self.inside) + ff.setNumber(tag, 'VOut', self.outside) + ff.setNumber(tag, 'XCenter', self.center[0]) + ff.setNumber(tag, 'YCenter', self.center[1]) + ff.setNumber(tag, 'ZCenter', 0) + return 'F{tag}' + + +@dataclass(frozen=True) +class Min(Field): + + f1: AsField + f2: AsField + + def getexpr(self, ff: FieldFactory) -> str: + s1 = as_field(self.f1).getexpr(ff) + s2 = as_field(self.f2).getexpr(ff) + return f'min({s1}, {s2})' + + +@dataclass(frozen=True) +class Max(Field): + + f1: AsField + f2: AsField + + def getexpr(self, ff: FieldFactory) -> str: + s1 = as_field(self.f1).getexpr(ff) + s2 = as_field(self.f2).getexpr(ff) + return f'max({s1}, {s2})' + + +def set_background(ff: FieldFactory, elemsize: AsField) -> None: + F = as_field(elemsize).getexpr(ff) + tag = ff.add('MathEval') + ff.setString(tag, 'F', F) + ff.setAsBackgroundMesh(tag) + + +# vim:sw=4:sts=4:et diff --git a/nutils/cs/_gmsh_stub.py b/nutils/cs/_gmsh_stub.py new file mode 100644 index 000000000..8dabac96d --- /dev/null +++ b/nutils/cs/_gmsh_stub.py @@ -0,0 +1,53 @@ +# Stub file for the used portions of the gmsh API + +from typing import Tuple, Protocol, NewType, NamedTuple, Sequence +from math import pi + +Tag = NewType('Tag', int) +DimTags = Sequence[Tuple[int, Tag]] + +class Affine(NamedTuple): + xx: float; xy: float; xz: float; dx: float + yx: float; yy: float; yz: float; dy: float + zx: float; zy: float; zz: float; dz: float + u1: float; u2: float; u3: float; u4: float + @classmethod + def shift(cls, dx: float, dy: float, dz: float): + return cls(1., 0., 0., dx, 0., 1., 0., dy, 0., 0., 1., dz, 0., 0., 0., 1.) + +class OCC(Protocol): + # https://gmsh.info/doc/texinfo/gmsh.html#Namespace-gmsh_002fmodel_002focc + def addRectangle(self, x: float, y: float, z: float, dx: float, dy: float, tag: Tag = Tag(-1), roundedRadius: float = 0.) -> Tag: ... + def addDisk(self, xc: float, yc: float, zc: float, rx: float, ry: float) -> Tag: ... + def addPoint(self, x: float, y: float, z: float, meshSize: float = 0., tag: Tag = Tag(-1)) -> Tag: ... + def addLine(self, startTag: Tag, endTag: Tag, tag: Tag = Tag(-1)) -> Tag: ... + def addCircleArc(self, startTag: Tag, centerTag: Tag, endTag: Tag, tag: Tag = Tag(-1), nx: float = 0., ny: float = 0., nz: float = 0.) -> Tag: ... + def addCurveLoop(self, curveTags: Sequence[Tag], tag: Tag = Tag(-1)) -> Tag: ... + def addPlaneSurface(self, wireTags: Sequence[Tag], tag: Tag = Tag(-1)) -> Tag: ... + def addBox(self, x: float, y: float, z: float, dx: float, dy: float, dz: float) -> Tag: ... + def addCylinder(self, x: float, y: float, z: float, dx: float, dy: float, dz: float, r: float) -> Tag: ... + def addSphere(self, xc: float, yc: float, zc: float, radius: float, tag: Tag = Tag(-1), angle1: float = pi/2, angle2: float = pi/2, angle3: float = 2*pi) -> Tag: ... + def addWire(self, curveTags: Sequence[Tag], tag: Tag = Tag(-1), checkClosed: bool = False) -> Tag: ... + def addPipe(self, dimTags: DimTags, wiretag: Tag, trihedron: str = "") -> DimTags: ... + def rotate(self, dimTags: DimTags, x: float, y: float, z: float, ax: float, ay: float, az: float, angle: float) -> Tag: ... + def fuse(self, objectDimTags: DimTags, toolDimTags: DimTags, tag: Tag = Tag(-1), removeObject: bool = True, removeTool: bool = True) -> Tuple[DimTags, Sequence[DimTags]]: ... + def cut(self, objectDimTags: DimTags, toolDimTags: DimTags, tag: Tag = Tag(-1), removeObject: bool = True, removeTool: bool = True) -> Tuple[DimTags, Sequence[DimTags]]: ... + def intersect(self, objectDimTags: DimTags, toolDimTags: DimTags, tag: Tag = Tag(-1), removeObject: bool = True, removeTool: bool = True) -> Tuple[DimTags, Sequence[DimTags]]: ... + def revolve(self, dimTags: DimTags, x: float, y: float, z: float, ax: float, ay: float, az: float, angle: float) -> DimTags: ... + def translate(self, dimTags: DimTags, dx: float, dy: float, dz: float) -> None: ... + def fragment(self, objectDimTags: DimTags, toolDimTags: Sequence[Tag], tag: Tag = Tag(-1), removeObject: bool = True, removeTool: bool = True) -> Tuple[Sequence[Tag],Sequence[DimTags]]: ... + def synchronize(self) -> None: ... + +class Mesh(Protocol): + # https://gmsh.info/doc/texinfo/gmsh.html#Namespace-gmsh_002fmodel_002fmesh + def setPeriodic(self, dim: int, tags: Sequence[Tag], tagsMaster: Sequence[Tag], affineTransform: 'Affine') -> None: ... + def generate(self, dim: int = 3) -> None: ... + +class Model(Protocol): + # https://gmsh.info/doc/texinfo/gmsh.html#Namespace-gmsh_002fmodel + occ: OCC + mesh: Mesh + def addPhysicalGroup(self, dim: int, tags: Sequence[Tag], tag: Tag = Tag(-1), name: str = '') -> Tag: ... + def setPhysicalName(self, dim: int, tag: Tag, name: str) -> None: ... + def getBoundary(self, dimTags: DimTags, combined: bool = True, oriented: bool = True, recursive: bool = False) -> DimTags: ... + def removeEntities(self, dimTags: DimTags, recursive: bool = False) -> None: ... diff --git a/nutils/cs/_shape.py b/nutils/cs/_shape.py new file mode 100644 index 000000000..7c866afad --- /dev/null +++ b/nutils/cs/_shape.py @@ -0,0 +1,485 @@ +from ._gmsh_stub import OCC, Mesh, Model, Tag, DimTags, Affine +from ._axes import Orientation +from ._util import overdimensioned +from typing import Tuple, Dict, Optional, Iterable, Sequence, List, Mapping, Set +from abc import ABC, abstractmethod +import numpy, re + + +Tags = Set[Tag] +XY = Tuple[float,float] +XYZ = Tuple[float,float,float] +Fragments = Dict['Shape', Tuple[Tags, Dict[str, Tags]]] +Entities = Mapping[str,'Entity'] + + +class Entity(ABC): + + def __init__(self, ndims: int): + self.ndims = ndims + + @abstractmethod + def get_shapes(self) -> Iterable['Shape']: + ... + + @abstractmethod + def select(self, fragments: Fragments) -> Tags: + ... + + +class Shape(Entity): + + def __init__(self, ndims: int, periodicity: Iterable[Tuple[str,str,Affine]] = ()): + self._periodicity = tuple(periodicity) + super().__init__(ndims) + + def get_shapes(self): + yield self + + def __sub__(self, other: 'Shape') -> 'BinaryOp': + return BinaryOp(self, other, 'cut') + + def __and__(self, other: 'Shape') -> 'BinaryOp': + return BinaryOp(self, other, 'intersect') + + def __or__(self, other: 'Shape') -> 'BinaryOp': + return BinaryOp(self, other, 'fuse') + + def extruded(self, segments: Sequence[Tuple[float,float,float]], **orientation_kwargs) -> 'Pipe': + '''Extruded 2D shape along 3D wire. + + The 2D `shape` is positioned in 3D space by translating to `origin` and + rotating in the directions of `xaxis` and `yaxis`. The shape is + subsequently extruded via a number of sections, each of which defines a + length, an x-curvature and a y-curvature. If both curvatures are zero then + the shape is linearly extruded over a distance of `length`. Otherwise, the + vector (xcurv, ycurv) defines the axis of rotation in the 2D plane, and its + length the curvature. Rotation follows the right-hand rule.''' + + orientation_kwargs.setdefault('origin', numpy.zeros(self.ndims+1)) + return Pipe(self, segments, Orientation(**orientation_kwargs)) + + def revolved(self, angle: float = 360., **orientation_kwargs) -> 'Revolved': + ''''Revolve 2D shape. + + The 2D `shape` is positioned in 3D space by translating to `origin` and + rotating in the directions of `xaxis` and `yaxis`. The shape is + subsequently rotated over its y-axis to form a (partially) revolved body. + In case the rotation `angle` is less than 360 degrees then boundary groups + 'front' and 'back' are added to the 2D shape's existing boundaries; + otherwise the revolved shape defines only the 2D boundary groups.''' + + orientation_kwargs.setdefault('origin', numpy.zeros(self.ndims+1)) + return Revolved(self, angle, Orientation(**orientation_kwargs)) + + def select(self, fragments: Fragments) -> Tags: + vtags, btags = fragments[self] + return vtags + + def make_periodic(self, mesh: Mesh, btags: Dict[str, Tags]) -> None: + for a, b, affinetrans in self._periodicity: + mesh.setPeriodic(self.ndims-1, sorted(btags[a]), sorted(btags[b]), affinetrans) + + @property + def boundary(self) -> 'Boundary': + return Boundary(self) + + @abstractmethod + def add_to(self, occ: OCC) -> DimTags: + ... + + @abstractmethod + def bnames(self, n: int) -> Iterable[str]: + ... + + +class Interval(Shape): + + def __init__(self, left: Optional[float] = None, right: Optional[float] = None, center: Optional[float] = None, length: Optional[float] = None, periodic: bool = False): + self.left, self.right, self.center, self.length = overdimensioned((left, 1, 0), (right, 0, 1), (center, .5, .5), (length, -1, 1)) + if self.length <= 0: + raise ValueError('negative interval') + self.periodic = periodic + super().__init__(ndims=1) + + def bnames(self, n): + assert n == 2 + return 'left', 'right' + + def add_to(self, occ: OCC) -> DimTags: + p1, p2 = [occ.addPoint(x, 0, 0) for x in (self.left, self.right)] + return (1, occ.addLine(p1, p2)), + + +class Point(Shape): + + def __init__(self, *p): + self.p = p + super().__init__(ndims=0) + + def bnames(self, n: int) -> Iterable[str]: + assert n == 0 + return () + + def add_to(self, occ: OCC) -> DimTags: + p = occ.addPoint(*self.p, *[0]*(3-len(self.p))) + return (0, p), + + +class Line(Shape): + + def __init__(self, p1, p2): + self.p1 = p1 + self.p2 = p2 + assert len(p1) == len(p2) + super().__init__(ndims=1) + + def bnames(self, n): + assert n == 2 + return 'left', 'right' + + def add_to(self, occ: OCC) -> DimTags: + p1, p2 = [occ.addPoint(*p, *[0]*(3-len(p))) for p in (self.p1, self.p2)] + return (self.ndims, occ.addLine(p1, p2)), + + +class Rectangle(Shape): + 'Rectangular domain' + + def __init__(self, x: Interval = Interval(0, 1), y: Interval = Interval(0, 1)): + self.x = x.left + self.dx = x.length + self.y = y.left + self.dy = y.length + periodicity = [(b, a, Affine.shift(*d*iv.length)) for d, (a, iv, b) in zip(numpy.eye(3), + [('left', x, 'right'), ('bottom', y, 'top')]) if iv.periodic] + super().__init__(ndims=2, periodicity=periodicity) + + def bnames(self, n): + assert n == 4 + return 'bottom', 'right', 'top', 'left' + + def add_to(self, occ: OCC) -> DimTags: + return (2, occ.addRectangle(x=self.x, y=self.y, z=0., dx=self.dx, dy=self.dy)), + + +class Circle(Shape): + 'Circular domain' + + def __init__(self, center: XY = (0., 0.), radius: float = 1.): + self.center = center + self.radius = radius + super().__init__(ndims=2) + + def bnames(self, n): + assert n == 1 + yield 'wall' + + def add_to(self, occ: OCC) -> DimTags: + return (2, occ.addDisk(*self.center, 0., rx=self.radius, ry=self.radius)), + + +class Ellipse(Shape): + 'Ellipsoidal domain' + + def __init__(self, center: XY = (0., 0.), width: float = 1., height: float = .5, angle: float = 0.): + self.center = center + self.width = width + self.height = height + self.angle = angle + super().__init__(ndims=2) + + def bnames(self, n): + assert n == 1 + yield 'wall' + + def add_to(self, occ: OCC) -> DimTags: + height, width, angle = (self.height, self.width, self.angle) if self.width > self.height \ + else (self.width, self.height, self.angle + 90) + tag = occ.addDisk(*self.center, 0., rx=width/2, ry=height/2) + occ.rotate([(2, tag)], *self.center, 0., ax=0., ay=0., az=1., angle=angle*numpy.pi/180) + return (2, tag), + + +class Path(Shape): + 'Arbitrarily shaped domain enclosed by straight and curved boundary segments' + + def __init__(self, vertices: Tuple[XY,...], curvatures: Optional[Tuple[float,...]] = None): + self.vertices = tuple(vertices) + assert all(len(v) == 2 for v in self.vertices) + self.curvatures = numpy.array(curvatures) if curvatures else numpy.zeros(len(self.vertices)) + assert len(self.curvatures) == len(self.vertices) + super().__init__(ndims=2) + + def bnames(self, n): + assert n == len(self.vertices) + return [f'segment{i}' for i in range(n)] + + def add_to(self, occ: OCC) -> DimTags: + points = [(v, occ.addPoint(*v, 0.)) for v in self.vertices] + points.append(points[0]) + lines = [occ.addLine(p1, p2) if k == 0 + else occ.addCircleArc(p1, occ.addPoint(*self._center(v1, v2, 1/k), 0.), p2) + for k, (v1, p1), (v2, p2) in zip(self.curvatures, points[:-1], points[1:])] + loop = occ.addCurveLoop(lines) + return (2, occ.addPlaneSurface([loop])), + + @staticmethod + def _center(v1: XY, v2: XY, r: float, nudge: float = 1e-7) -> XY: + cx, cy = numpy.add(v1, v2) / 2 + dx, dy = numpy.subtract(v1, v2) / 2 + r2 = dx**2 + dy**2 + D2 = r**2 / r2 - 1 + nudge + assert D2 > 0, f'invalid arc radius: {r} < {numpy.sqrt(r2)}' + D = numpy.copysign(numpy.sqrt(D2), r) + return cx + dy * D, cy - dx * D, + + +class Box(Shape): + 'Box' + + def __init__(self, x: Interval = Interval(0, 1), y: Interval = Interval(0, 1), z: Interval = Interval(0, 1)): + self.x = x.left + self.dx = x.length + self.y = y.left + self.dy = y.length + self.z = z.left + self.dz = z.length + periodicity = [(b, a, Affine.shift(*d*iv.length)) for d, (a, iv, b) in zip(numpy.eye(3), + [('left', x, 'right'), ('bottom', y, 'top'), ('front', z, 'back')]) if iv.periodic] + super().__init__(ndims=3, periodicity=periodicity) + + def bnames(self, n): + assert n == 6 + return 'left', 'right', 'bottom', 'top', 'front', 'back' + + def add_to(self, occ: OCC) -> DimTags: + return (3, occ.addBox(x=self.x, y=self.y, z=self.z, dx=self.dx, dy=self.dy, dz=self.dz)), + + +class Sphere(Shape): + 'Sphere' + + def __init__(self, center: XYZ = (0., 0., 0., ), radius: float = 1.): + self.center = center + self.radius = radius + super().__init__(ndims=3) + + def bnames(self, n): + assert n == 1 + return 'wall', + + def add_to(self, occ: OCC) -> DimTags: + return (3, occ.addSphere(*self.center, self.radius)), + + +class Cylinder(Shape): + 'Cylinder' + + def __init__(self, front: XYZ = (0.,0.,0.), back: XYZ = (0.,0.,1.), radius: float = 1., periodic: bool = False): + self.center = front + self.axis = back[0] - front[0], back[1] - front[1], back[2] - front[2] + self.radius = radius + super().__init__(ndims=3, periodicity=[('back', 'front', Affine.shift(*self.axis))] if periodic else ()) + + def bnames(self, n): + assert n == 3 + return 'side', 'back', 'front' + + def add_to(self, occ: OCC) -> DimTags: + return (3, occ.addCylinder(*self.center, *self.axis, self.radius)), + + +class BinaryOp(Shape): + + def __init__(self, shape1: Shape, shape2: Shape, op: str): + self.shape1 = shape1 + self.shape2 = shape2 + self.op = op + assert shape2.ndims == shape1.ndims + super().__init__(shape1.ndims) + + def bnames(self, n): + return [f'section{i}' for i in range(n)] + + def add_to(self, occ: OCC) -> DimTags: + op = getattr(occ, self.op) + return op(objectDimTags=self.shape1.add_to(occ), toolDimTags=self.shape2.add_to(occ))[0] + + +class Revolved(Shape): + + def __init__(self, shape: Shape, angle: float, orientation: Orientation): + assert orientation.ndims == shape.ndims + 1 + self.shape = shape + self.front = orientation + self.angle = float(angle) * numpy.pi / 180 + super().__init__(ndims=orientation.ndims) + + def bnames(self, n): + partial = self.angle < 2 * numpy.pi + if partial: + yield 'back' + n -= 2 + for bname in self.shape.bnames(n): + yield f'side-{bname}' + if partial: + yield 'front' + + def add_to(self, occ: OCC) -> DimTags: + front = self.shape.add_to(occ) + self.front.orient(occ, front) + axes, origin = self.front.axes.as_3(self.front.origin) + iaxis = {2: 2, 3: 1}[self.ndims] # TODO: allow variation of revolution axis in 3D + dimtags = occ.revolve(front, *origin, *axes[iaxis], -self.angle) + return [(dim, tag) for dim, tag in dimtags if dim == self.ndims] + + +class Pipe(Shape): + + def __init__(self, shape: Shape, segments: Sequence[Tuple[float,float,float]], orientation: Orientation): + assert orientation.ndims == shape.ndims + 1 + self.shape = shape + self.nsegments = len(segments) + + self.front = orientation + vertices = [orientation.origin] + midpoints: List[Optional[XYZ]] = [] + for length, *curvature in segments: + if not any(curvature): + midpoints.append(None) + orientation = Orientation(orientation.origin + length * orientation.axes[-1], axes=orientation.axes) + else: + if shape.ndims == 1: + kx, = curvature + radius = numpy.array([-1/kx]) + rotation = kx * length + elif shape.ndims == 2: + kx, ky = curvature + radius = numpy.divide([ky, -kx], kx**2 + ky**2) + rotation = numpy.multiply(curvature, length) @ orientation.axes[:-1] + else: + raise NotImplementedError + center = orientation.origin + radius @ orientation.axes[:-1] + midpoints.append(center) + axes = orientation.axes.rotate(rotation) + orientation = Orientation(center - radius @ axes[:-1], axes=axes) + vertices.append(orientation.origin) + self.midpoints = tuple(midpoints) + self.vertices = tuple(vertices) + self.back = orientation + + super().__init__(ndims=orientation.ndims) + + def bnames(self, n): + nb, rem = divmod(n-2, self.nsegments) + assert not rem + if self.ndims == 2: + bnames = [f'segment{i}-{bname}' for i in range(self.nsegments) for bname in self.shape.bnames(nb)] + bnames.insert(1, 'front') + bnames.append('back') + elif self.ndims == 3: + bnames = [f'segment{i}-{bname}' for bname in self.shape.bnames(nb) for i in range(self.nsegments)] + bnames.insert(0, 'front') + bnames.append('back') + else: + raise NotImplementedError + return bnames + + def add_to(self, occ: OCC) -> DimTags: + z = (0,) * (3 - self.ndims) + points = [occ.addPoint(*v, *z) for v in self.vertices] + segments = [occ.addLine(p1, p2) if v is None + else occ.addCircleArc(p1, occ.addPoint(*v, *z), p2) for p1, v, p2 in zip(points, self.midpoints, points[1:])] + wire_tag = occ.addWire(segments) + front = self.shape.add_to(occ) + self.front.orient(occ, front) + return occ.addPipe(front, wire_tag) + + +class Boundary(Entity): + + def __init__(self, parent: Shape, patterns = ()): + self.parent = parent + self.patterns = patterns + super().__init__(parent.ndims - 1) + + def __getitem__(self, item: str) -> 'Boundary': + return Boundary(self.parent, (*self.patterns, re.compile(item))) + + def get_shapes(self): + return self.parent.get_shapes() + + def select(self, fragments: Fragments) -> Tags: + vtags, btags = fragments[self.parent] + s = set.union(*[items for bname, items in btags.items() if all(pattern.fullmatch(bname) for pattern in self.patterns)]) + if not s: + raise ValueError(f'{self.parent} does not have a boundary {", ".join(p.pattern for p in self.patterns)}') + return s + + +class Skeleton(Entity): + + def get_shapes(self): + return () + + def select(self, fragments: Fragments) -> Tags: + return set.union(*[items for vtags, btags in fragments.values() for items in btags.values()]) + + +def generate_mesh(model: Model, mapping: Entities) -> None: + + shapes = tuple(dict.fromkeys(shape for entity in mapping.values() for shape in entity.get_shapes())) # stable unique via dict + shape_tags = [shape.add_to(model.occ) for shape in shapes] # create all shapes before sync + + model.occ.synchronize() # required for getBoundary + + objectDimTags: List[Tuple[int, Tag]] = [] + slices = [] + a = 0 + for dimtags in shape_tags: + objectDimTags.extend(dimtags) + b = len(objectDimTags) + vslice = slice(a, b) + objectDimTags.extend(model.getBoundary(dimtags, oriented=False)) + a = len(objectDimTags) + bslice = slice(b, a) + slices.append((vslice, bslice)) + _, fragment_map = model.occ.fragment(objectDimTags=objectDimTags, toolDimTags=[], removeObject=False) + assert len(fragment_map) == a + + model.occ.synchronize() + + # setting fragment's removeObject=True has a tendency to remove (boundary) + # entities that are still in use, so we remove unused entities manually + # instead + remove = set(objectDimTags) + for dimtags in fragment_map: + remove.difference_update(dimtags) + if remove: + model.removeEntities(sorted(remove)) + + fragments = {} + for shape, (vslice, bslice) in zip(shapes, slices): + + vtags = _tags([dimtag for dimtags in fragment_map[vslice] for dimtag in dimtags], shape.ndims) + btagslist = [_tags(dimtags, shape.ndims-1) for dimtags in fragment_map[bslice]] + + btags = dict(zip(shape.bnames(len(btagslist)), btagslist)) + shape.make_periodic(model.mesh, btags) + + fragments[shape] = vtags, btags + + for name, shape in mapping.items(): + tag = model.addPhysicalGroup(shape.ndims, sorted(shape.select(fragments))) + model.setPhysicalName(dim=shape.ndims, tag=tag, name=name) + + model.mesh.generate(max(shape.ndims for shape in shapes)) + + +def _tags(dimtags: DimTags, expect_dim: int) -> Tags: + assert all(dim == expect_dim for dim, tag in dimtags) + return {tag for dim, tag in dimtags} + + +# vim:sw=4:sts=4:et diff --git a/nutils/cs/_util.py b/nutils/cs/_util.py new file mode 100644 index 000000000..9cdb2ab5a --- /dev/null +++ b/nutils/cs/_util.py @@ -0,0 +1,15 @@ +import numpy + +def overdimensioned(*args): + ncoeffs = len(args[0]) - 1 + matrix = [] + values = [] + for v, *coeffs in args: + assert len(coeffs) == ncoeffs + if v is not None: + matrix.append(coeffs) + values.append(v) + if len(values) != ncoeffs: + raise ValueError(f'exactly {ncoeffs} arguments should be specified') + mat = numpy.linalg.solve(matrix, values).T + return [mat @ coeffs for v, *coeffs in args] diff --git a/nutils/mesh.py b/nutils/mesh.py index 3e9bc0f2d..b215d204c 100644 --- a/nutils/mesh.py +++ b/nutils/mesh.py @@ -6,7 +6,7 @@ provided at this point. """ -from . import topology, function, _util as util, element, numeric, transform, transformseq, warnings, types, cache +from . import topology, function, _util as util, element, numeric, transform, transformseq, warnings, types, cache, cs from .elementseq import References from .transform import TransformItem from .topology import Topology @@ -474,6 +474,22 @@ def gmsh(fname, name='gmsh', *, space='X'): return simplex(name=name, **parsegmsh(f), space=space) +def parsecsgmsh(entities: cs.Entities, elemsize: cs.AsField, order: int = 1): + from tempfile import mkstemp + fid, fname = mkstemp(suffix='.msh') + try: + os.close(fid) # release file for writing by gmsh (windows) + cs.write(fname, entities, elemsize, order) + with open(fname, 'rb') as f: + return parsegmsh(f) + finally: + os.unlink(fname) + + +def csgmsh(entities: cs.Entities, elemsize: cs.AsField, order: int = 1, name='csgmsh', space='X'): + return simplex(name=name, space=space, **parsecsgmsh(entities, elemsize, order)) + + def simplex(nodes, cnodes, coords, tags, btags, ptags, name='simplex', *, space='X'): '''Simplex topology. diff --git a/tests/test_cs.py b/tests/test_cs.py new file mode 100644 index 000000000..d876f09da --- /dev/null +++ b/tests/test_cs.py @@ -0,0 +1,406 @@ +from nutils import mesh, function, cs +import unittest +import numpy + + +class TestShape(unittest.TestCase): + + def assertAreaCentroid(self, topo, geom, *, mass, centroid, places, degree=1): + J = function.J(geom) + _mass, _moment = topo.sample('gauss', 2).integrate([J, geom * J], degree=degree*2) + self.assertAlmostEqual(_mass, mass, places=places) + numpy.testing.assert_almost_equal(actual=_moment/mass, desired=centroid, decimal=places) + + def test_rectangle(self): + 'Square with 1 Date: Thu, 24 Aug 2023 16:01:23 +0200 Subject: [PATCH 2/5] add csgmsh topo to platewithhole example --- examples/platewithhole.py | 56 +++++++++++++++++++++++++++++++-------- 1 file changed, 45 insertions(+), 11 deletions(-) diff --git a/examples/platewithhole.py b/examples/platewithhole.py index 5f8845858..8bd6904c8 100644 --- a/examples/platewithhole.py +++ b/examples/platewithhole.py @@ -1,4 +1,4 @@ -from nutils import mesh, function, export, testing +from nutils import mesh, function, export, testing, cs from nutils.solver import System from nutils.expression_v2 import Namespace from dataclasses import dataclass @@ -43,8 +43,8 @@ def generate(self, radius): @dataclass -class NURBS: - '''Non-Uniform Radional B-Splines +class IGA: + '''Isogeometric Analysis Generate a 1x2 structured topology, map it using quadratic NURBS to a square domain with circular cut-out, and refine several times before @@ -81,7 +81,33 @@ def generate(self, radius): return topo.withboundary(hole='left', sym='top,bottom', far='right'), geom, nurbsbasis, 5 -def main(mode: Union[FCM, NURBS] = NURBS(), +@dataclass +class TRI: + '''Triangulation + + Generate a triangulated mesh with second order geometry using gmsh. + + Parameters + ---------- + elemsize + Target element size. + degree + Polynomial degree for the basis. + ''' + + elemsize: int = .1 + degree: int = 2 + + def generate(self, radius): + rect = cs.Rectangle() + circ = cs.Circle(radius=radius) + shapes = dict(dom=rect-circ, sym=rect.boundary['left|bottom'], far=rect.boundary['right|top'], hole=circ.boundary) + topo, geom = mesh.csgmsh(shapes, elemsize=self.elemsize, order=2) + basis = topo.basis('std', self.degree) + return topo, geom, basis, self.degree + + +def main(mode: Union[FCM, IGA, TRI] = IGA(), radius: float = .5, traction: float = .1, poisson: float = .3): @@ -96,7 +122,7 @@ def main(mode: Union[FCM, NURBS] = NURBS(), The script can be run in two modes: by specifying `mode=FCM`, the circular hole is cut out of a regular finite element mesh by means of the Finite Cell - Method; by specifying `mode=NURBS` a Non-Uniform Rational BSpline geometry is + Method; by specifying `mode=IGA` a Non-Uniform Rational BSpline geometry is created to map a regular domain onto the desired shape. Either mode supports sub-parameters which can be specified from the command-line by attaching them in curly braces (e.g. `FCM{nelems=20,degree=1}`). @@ -104,7 +130,8 @@ def main(mode: Union[FCM, NURBS] = NURBS(), Parameters ---------- mode - Discretization strategy: FCM (Finite Cell Method) or NURBS. + Discretization strategy: FCM (Finite Cell Method), IGA (Isogeometric + Analysis) or TRI (triangulation). radius Cut-out radius. traction @@ -137,10 +164,10 @@ def main(mode: Union[FCM, NURBS] = NURBS(), log.info('hole radius exact up to L2 error {:.2e}'.format(radiuserr)) sqr = topo.boundary['sym'].integral('(u_i n_i)^2 dS' @ ns, degree=degree*2) - cons = System(sqr, trial='u').solve_constraints(droptol=1e-15) + cons = System(sqr, trial='u').solve_constraints(droptol=1e-10) sqr = topo.boundary['far'].integral('du_k du_k dS' @ ns, degree=20) - cons = System(sqr, trial='u').solve_constraints(droptol=1e-15, constrain=cons) + cons = System(sqr, trial='u').solve_constraints(droptol=1e-10, constrain=cons) res = topo.integral('∇_j(v_i) σ_ij dV' @ ns, degree=degree*2) args = System(res, trial='u', test='v').solve(constrain=cons) @@ -193,7 +220,7 @@ def test_mixed(self): fyhZkYI=''') def test_nurbs0(self): - err, cons, args = main(mode=NURBS(nrefine=0)) + err, cons, args = main(mode=IGA(nrefine=0)) with self.subTest('l2-error'): self.assertAlmostEqual(err[0], .00200, places=5) with self.subTest('h1-error'): @@ -206,7 +233,7 @@ def test_nurbs0(self): eNpjYJh07qLhhnOTjb0vTDdmAAKVcy/1u85lGYforQDzFc6pGSedlzd+eP4ykA8AvkQRaA==''') def test_nurbs2(self): - err, cons, args = main(mode=NURBS(nrefine=2)) + err, cons, args = main(mode=IGA(nrefine=2)) with self.subTest('l2-error'): self.assertAlmostEqual(err[0], .00009, places=5) with self.subTest('h1-error'): @@ -222,10 +249,17 @@ def test_nurbs2(self): n2c23nZe3djqQqpx88XNxrOv7gOr0zwXZeBxztro/bmnRp7nVY1zgTjvvIXxSaBfnl3YYbzmygmgOgDU Imlr''') + def test_tri(self): + err, cons, args = main(mode=TRI(elemsize=.5)) + with self.subTest('l2-error'): + self.assertAlmostEqual(err[0], .00053, places=5) + with self.subTest('h1-error'): + self.assertAlmostEqual(err[1], .0131, places=4) + if __name__ == '__main__': from nutils import cli cli.run(main) -# example:tags=elasticity,FCM,NURBS +# example:tags=elasticity,FCM,IGA From 79bba97e6c213bf1891d9c80701c37a250d29b27 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Thu, 24 Aug 2023 16:24:09 +0200 Subject: [PATCH 3/5] add csgmsh, multipatch topo to adaptivity example --- examples/adaptivity.py | 82 +++++++++++++++++++++++++++--------------- 1 file changed, 54 insertions(+), 28 deletions(-) diff --git a/examples/adaptivity.py b/examples/adaptivity.py index 71d468e1c..643eb6ff9 100644 --- a/examples/adaptivity.py +++ b/examples/adaptivity.py @@ -1,11 +1,44 @@ -from nutils import mesh, function, export, testing +from nutils import mesh, function, export, testing, cs from nutils.solver import System from nutils.expression_v2 import Namespace +from dataclasses import dataclass +from typing import Union import numpy import treelog -def main(etype: str = 'square', +@dataclass +class structured: + + def generate(self): + topo, geom = mesh.rectilinear([2,2]) + quadrant = topo[1:,:1].withboundary(corner='left,top') + return topo - quadrant, (geom - 1) / 2 + + +@dataclass +class unstructured: + + elemsize: float = .15 + + def generate(self): + square = cs.Rectangle(cs.Interval(-.5, .5), cs.Interval(-.5, .5)) + quadrant = cs.Rectangle(cs.Interval(0, .5), cs.Interval(-.5, 0)) + shapes = dict(dom=square-quadrant, corner=quadrant.boundary) + return mesh.csgmsh(shapes, elemsize=self.elemsize) + + +@dataclass +class multipatch: + + def generate(self): + vertices = [0,-.5], [-.5,-.5], [0,0], [-.5,.5], [.5,0], [.5,.5] + patches = [0,1,2,3], [2,3,4,5] + topo, geom = mesh.multipatch(patchverts=vertices, patches=patches, nelems=1) + return topo.withboundary(corner='patch0-bottom,patch1-bottom'), geom + + +def main(basemesh: Union[structured, unstructured, multipatch] = structured(), btype: str = 'h-std', degree: int = 2, nrefine: int = 5): @@ -29,8 +62,8 @@ def main(etype: str = 'square', Parameters ---------- - etype - Type of elements (square/triangle/mixed). + basemesh + Initial mesh: structured, unstructured, or multipatch. btype Type of basis function (h/th-std/spline), with availability depending on the configured element type. @@ -38,15 +71,13 @@ def main(etype: str = 'square', Polynomial degree. nrefine Number of refinement steps to perform. + unstructured + Generate triangulated domain if True, or structured domain if False. ''' - domain, geom = mesh.unitsquare(2, etype) - geom -= .5 # shift domain center to origin - + domain, geom = basemesh.generate() x, y = geom exact = (x**2 + y**2)**(1/3) * numpy.cos(numpy.arctan2(y+x, y-x) * (2/3)) - selection = domain.select(exact, ischeme='gauss1') - domain = domain.subset(selection, newboundary='corner') linreg = LinearRegressor(bias=1) for irefine in treelog.iter.fraction('level', range(nrefine+1)): @@ -136,8 +167,8 @@ def offset(self): class test(testing.TestCase): - def test_square_quadratic(self): - error, u = main(nrefine=2) + def test_structured(self): + error, u = main(nrefine=2, basemesh=structured()) with self.subTest('degrees of freedom'): self.assertEqual(len(u), 149) with self.subTest('L2-error'): @@ -152,32 +183,27 @@ def test_square_quadratic(self): 7AYLvMPpsqGkCTPumzWf+qV92kKevjK36ozDP/FSnh1iteWiqWuf+oMaKuyKaC1i52rKPokiF2WLA/20 bya+ZCPbWKRPpvgFaedebw==''') - def test_triangle_quadratic(self): - error, u = main(nrefine=2, etype='triangle') + def test_unstructured(self): + error, u = main(nrefine=2, basemesh=unstructured(elemsize=.5)) with self.subTest('degrees of freedom'): - self.assertEqual(len(u), 98) + self.assertEqual(len(u), 97) with self.subTest('L2-error'): - self.assertAlmostEqual(error[0], 0.00138, places=5) + self.assertAlmostEqual(error[0], 0.00095, places=5) with self.subTest('H1-error'): - self.assertAlmostEqual(error[1], 0.05326, places=5) - with self.subTest('left-hand side'): - self.assertAlmostEqual64(u, ''' - eNprMV1oesqU2VTO1Nbko6myWbhpq+kckwST90avjRgYzptYm+YYMwBBk3GQWavZb1NXs2+mm83um1WY - bQbyXYEiQWbKZjNM7wJVzjBlYICoPW8CMiXH+LXRR9NwoPkg82xN5IB2MZu2mGabSBnnAbGscYEJj3GV - YQAQg/TVGfaA7RI0BsErRjeNeowDgDQPmF9gkmciaJxtArGjzrAKCGWNpYAQAL0kOBE=''') + self.assertAlmostEqual(error[1], 0.04006, places=5) - def test_mixed_linear(self): - error, u = main(nrefine=2, etype='mixed', degree=1) + def test_multipatch(self): + error, u = main(nrefine=2, basemesh=multipatch()) with self.subTest('degrees of freedom'): - self.assertEqual(len(u), 34) + self.assertEqual(len(u), 93) with self.subTest('L2-error'): - self.assertAlmostEqual(error[0], 0.00450, places=5) + self.assertAlmostEqual(error[0], 0.00128, places=5) with self.subTest('H1-error'): - self.assertAlmostEqual(error[1], 0.11692, places=5) + self.assertAlmostEqual(error[1], 0.05662, places=5) with self.subTest('left-hand side'): self.assertAlmostEqual64(u, ''' - eNprMT1u6mQyxUTRzMCUAQhazL6b3jNrMYPxp5iA5FtMD+lcMgDxHa4aXzS+6HDV+fKO85cMnC8zMBzS - AQDBThbY''') + eNpjYHhnZGfCZNpqysAgYJxuYmK6C8gqM35ksspU2QyZhZBF6DhnUGuUZPzJON5ktclbEyWgiInRVuOJ + QNlE06WmD02FgPpajaabuJpuMuUy8zZrMdsMFCk24TKda/rN1MWsDsi/j1UNpjnodgEAAk420A==''') if __name__ == '__main__': From de3920aef9c6f0437bf5be043999383ecf8e1342 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Fri, 25 Aug 2023 16:18:39 +0200 Subject: [PATCH 4/5] add "cs" dependency --- .github/workflows/test.yaml | 10 ++++++++-- pyproject.toml | 1 + 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml index 937ddc46e..d25701499 100644 --- a/.github/workflows/test.yaml +++ b/.github/workflows/test.yaml @@ -89,6 +89,9 @@ jobs: - name: Install Graphviz if: ${{ matrix.os == 'ubuntu-latest' }} run: sudo apt install -y graphviz + - name: Install LibGLU + if: ${{ matrix.os == 'ubuntu-latest' }} + run: sudo apt install -y libglu1 - name: Install Nutils and dependencies id: install env: @@ -97,7 +100,7 @@ jobs: python -um pip install --upgrade --upgrade-strategy eager wheel python -um pip install --upgrade --upgrade-strategy eager numpy$_numpy_version # Install Nutils from `dist` dir created in job `build-python-package`. - python -um pip install "$_wheel[import_gmsh,export_mpl]" + python -um pip install "$_wheel[import_gmsh,export_mpl,cs]" - name: Install Scipy if: ${{ matrix.matrix-backend == 'scipy' }} run: python -um pip install --upgrade scipy @@ -171,12 +174,15 @@ jobs: with: name: python-package path: dist/ + - name: Install LibGLU + if: ${{ matrix.os == 'ubuntu' }} + run: sudo apt install -y libglu1 - name: Install Nutils and dependencies id: install run: | python -um pip install --upgrade --upgrade-strategy eager wheel # Install Nutils from `dist` dir created in job `build-python-package`. - python -um pip install "$_wheel[matrix_scipy,export_mpl]" + python -um pip install "$_wheel[matrix_scipy,export_mpl,cs]" - name: Test run: python -um unittest discover -b -q -t . -s examples test-sphinx: diff --git a/pyproject.toml b/pyproject.toml index ab7a4ed35..fb6e1f0e8 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -27,6 +27,7 @@ export_mpl = ["matplotlib >=3.3,<4"] matrix_mkl = ["mkl"] matrix_scipy = ["scipy >=0.13,<2"] import_gmsh = ["meshio >=4,<6"] +cs = ["gmsh>=4", "meshio>=5"] [build-system] requires = ["flit_core >=3.2,<4"] From 4f0aa08ef9c39bea200309c2c956884a1f758e5b Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Fri, 12 Jul 2024 10:27:29 +0200 Subject: [PATCH 5/5] wip --- nutils/cs/__init__.py | 6 +- nutils/cs/_field.py | 183 ++++++++++++++++++++++++++++-------------- nutils/cs/_shape.py | 11 ++- 3 files changed, 132 insertions(+), 68 deletions(-) diff --git a/nutils/cs/__init__.py b/nutils/cs/__init__.py index f4b156967..f0820bbae 100644 --- a/nutils/cs/__init__.py +++ b/nutils/cs/__init__.py @@ -22,6 +22,8 @@ AsField, Ball, LocalRefinement, + Distance, + Threshold, Max, Min, set_background, @@ -47,8 +49,8 @@ def write(fname: str, entities: Entities, elemsize: AsField, order: int = 1) -> gmsh.logger.start() try: - set_background(gmsh.model.mesh.field, elemsize) - generate_mesh(gmsh.model, entities) + #set_background(gmsh.model.mesh.field, elemsize) + generate_mesh(gmsh.model, entities, elemsize) gmsh.write(fname) finally: with treelog.context('gmsh'): diff --git a/nutils/cs/_field.py b/nutils/cs/_field.py index 07392a3e9..8ca998c30 100644 --- a/nutils/cs/_field.py +++ b/nutils/cs/_field.py @@ -1,4 +1,5 @@ from typing import Tuple, Union, Protocol, runtime_checkable, NewType +import math from dataclasses import dataclass @@ -26,80 +27,140 @@ def setAsBackgroundMesh(self, tag: Tag) -> None: class Field(Protocol): 'A scalar field over a coordinate system.' - def getexpr(self, ff: FieldFactory) -> str: + def gettag(self, ff: FieldFactory, fragments) -> str: ... def __add__(self, other: AsField) -> 'Field': - return FieldOp(self, as_field(other), '+') + return MathEval('({})+({})', self, as_field(other)) def __radd__(self, other: AsField) -> 'Field': - return FieldOp(as_field(other), self, '+') + return MathEval('({})+({})', as_field(other), self) def __sub__(self, other: AsField) -> 'Field': - return FieldOp(self, as_field(other), '-') + return MathEval('({})-({})', self, as_field(other)) def __rsub__(self, other: AsField) -> 'Field': - return FieldOp(as_field(other), self, '-') + return MathEval('({})-({})', as_field(other), self) def __mul__(self, other: AsField) -> 'Field': - return FieldOp(self, as_field(other), '*') + return MathEval('({})*({})', as_field(other), self) def __rmul__(self, other: AsField) -> 'Field': - return FieldOp(as_field(other), self, '*') + return MathEval('({})*({})', self, as_field(other)) def __truediv__(self, other: AsField) -> 'Field': - return FieldOp(self, as_field(other), '/') + return MathEval('({})/({})', as_field(other), self) + def __pow__(self, other: AsField) -> 'Field': + return MathEval('({})^({})', as_field(other), self) def as_field(f: AsField) -> Field: if isinstance(f, Field): return f if isinstance(f, (int, float)): - return Constant(f) + return MathEval(str(f)) raise ValueError(f'cannot interpret {f!r} as a field') -@dataclass(frozen=True) -class FieldOp(Field): - - a: Field - b: Field - op: str +class MathEval(Field): - def getexpr(self, ff: FieldFactory) -> str: - a = self.a.getexpr(ff) - b = self.b.getexpr(ff) - return f'({a}){self.op}({b})' + def __init__(self, fmt, *args): + self.fmt = fmt + self.args = args + def getexpr(self, ff: FieldFactory, fragments) -> str: + return self.fmt.format(*[arg.getexpr(ff, fragments) if isinstance(arg, MathEval) else f'F{arg.gettag(ff, fragments)}' for arg in self.args]) -@dataclass(frozen=True) -class Constant(Field): - 'Constant element size' + def gettag(self, ff: FieldFactory, fragments): + expr = self.getexpr(ff, fragments) + tag = ff.add('MathEval') + print(tag, '->', expr) + ff.setString(tag, 'F', expr) + return tag - size: float - def getexpr(self, ff: FieldFactory) -> str: - return str(self.size) +x = MathEval('x') +y = MathEval('y') +z = MathEval('z') @dataclass(frozen=True) -class Coord(Field): +class LocalRefinement(Field): + 'Refine elements according to a gaussian bell.' - coord: str + distance: Field + radius: float + factor: float = 2. - def getexpr(self, ff: FieldFactory) -> str: - return self.coord + # d = 0 -> X + # d = radius -> Y + # d = inf -> 1 -x = Coord('x') -y = Coord('y') -z = Coord('z') + # scale = 1 - (1-X) * ((1-X)/(1-Y))^(-d^2/radius^2) + # X = 1/factor + # Y = 2/(factor+1) -@dataclass(frozen=True) -class LocalRefinement(Field): - 'Refine elements according to a gaussian bell.' + # (1-X)/(1-Y) = (1-1/factor)/(1-2/(factor+1)) = (factor+1)/factor = 1+1/factor + + # scale = 1 - (1-1/factor) * (1+1/factor)^(-d^2/radius^2) + + def gettag(self, ff: FieldFactory, fragments) -> str: + c = 1 / self.factor + d = self.distance.gettag(ff, fragments) + expr = f'1-{1-c}*{1+c}^(-(({d})/{self.radius})^2)' + #expr = f'1-.5*exp(-({d})^2)' + #s = 1 - 1/self.factor + #expr = f'1-{s}*exp(-(({d})/({self.radius}))^2)' + print(expr) + return expr - center: Tuple[float,float] = (0., 0.) - radius: float = 1. - factor: float = 2. - def getexpr(self, ff: FieldFactory) -> str: - x, y = self.center - s = 1 - 1/self.factor - return f'1-{s}*exp(-((x-{x})^2+(y-{y})^2)/({self.radius})^2)' +class Distance(Field): + 'Refine elements uniformly inside a circle' + + def __init__(self, *entities, sampling: int = 20): + self.entities = entities + self.sampling = sampling + + def gettag(self, ff: FieldFactory, fragments) -> str: + tag = ff.add('Distance') + ff.setNumber(tag, 'NumPointsPerCurve', round(self.sampling)) # gmsh 34a4d3c613 + #ff.setNumber(tag, 'Sampling', round(self.sampling)) + surfaces = set() + curves = set() + points = set() + for entity in self.entities: + tags = entity.select(fragments) + if entity.ndims == 2: + surfaces.update(tags) + elif entity.ndims == 1: + curves.update(tags) + elif entity.ndims == 0: + points.update(tags) + else: + bla + if surfaces: + ff.setNumbers(tag, 'SurfacesList', sorted(surfaces)) + if curves: + ff.setNumbers(tag, 'CurvesList', sorted(curves)) + if points: + ff.setNumbers(tag, 'PointsList', sorted(points)) + return tag + + +@dataclass +class Threshold(Field): + + d: Field + dmin: float + dmax: float + vmin: float + vmax: float + sigmoid: bool = False + + def gettag(self, ff: FieldFactory, fragments): + tag = ff.add('Threshold') + ff.setNumber(tag, 'InField', self.d.gettag(ff, fragments)) + ff.setNumber(tag, 'DistMin', self.dmin) + ff.setNumber(tag, 'DistMax', self.dmax) + ff.setNumber(tag, 'SizeMin', self.vmin) + ff.setNumber(tag, 'SizeMax', self.vmax) + ff.setNumber(tag, 'Sigmoid', self.sigmoid) + return tag @dataclass(frozen=True) @@ -112,7 +173,7 @@ class Ball(Field): outside: float thickness: float = 0. - def getexpr(self, ff: FieldFactory) -> str: + def gettag(self, ff: FieldFactory, fragments) -> str: tag = ff.add('Ball') ff.setNumber(tag, 'Radius', self.radius) ff.setNumber(tag, 'Thickness', self.thickness) @@ -121,37 +182,35 @@ def getexpr(self, ff: FieldFactory) -> str: ff.setNumber(tag, 'XCenter', self.center[0]) ff.setNumber(tag, 'YCenter', self.center[1]) ff.setNumber(tag, 'ZCenter', 0) - return 'F{tag}' + return tag -@dataclass(frozen=True) class Min(Field): - f1: AsField - f2: AsField + def __init__(self, *fields): + self.fields = fields - def getexpr(self, ff: FieldFactory) -> str: - s1 = as_field(self.f1).getexpr(ff) - s2 = as_field(self.f2).getexpr(ff) - return f'min({s1}, {s2})' + def gettag(self, ff: FieldFactory, fragments) -> str: + tags = [as_field(f).gettag(ff, fragments) for f in self.fields] + tag = ff.add('Min') + ff.setNumbers(tag, 'FieldsList', tags) + return tag -@dataclass(frozen=True) class Max(Field): - f1: AsField - f2: AsField + def __init__(self, *fields): + self.fields = fields - def getexpr(self, ff: FieldFactory) -> str: - s1 = as_field(self.f1).getexpr(ff) - s2 = as_field(self.f2).getexpr(ff) - return f'max({s1}, {s2})' + def gettag(self, ff: FieldFactory, fragments) -> str: + tags = [as_field(f).gettag(ff, fragments) for f in self.fields] + tag = ff.add('Max') + ff.setNumbers(tag, 'FieldsList', tags) + return tag -def set_background(ff: FieldFactory, elemsize: AsField) -> None: - F = as_field(elemsize).getexpr(ff) - tag = ff.add('MathEval') - ff.setString(tag, 'F', F) +def set_background(ff: FieldFactory, elemsize: AsField, fragments) -> None: + tag = as_field(elemsize).gettag(ff, fragments) ff.setAsBackgroundMesh(tag) diff --git a/nutils/cs/_shape.py b/nutils/cs/_shape.py index 7c866afad..2eb8b5b72 100644 --- a/nutils/cs/_shape.py +++ b/nutils/cs/_shape.py @@ -1,6 +1,7 @@ from ._gmsh_stub import OCC, Mesh, Model, Tag, DimTags, Affine from ._axes import Orientation from ._util import overdimensioned +from ._field import set_background from typing import Tuple, Dict, Optional, Iterable, Sequence, List, Mapping, Set from abc import ABC, abstractmethod import numpy, re @@ -427,7 +428,7 @@ def select(self, fragments: Fragments) -> Tags: return set.union(*[items for vtags, btags in fragments.values() for items in btags.values()]) -def generate_mesh(model: Model, mapping: Entities) -> None: +def generate_mesh(model: Model, mapping: Entities, elemsize) -> None: shapes = tuple(dict.fromkeys(shape for entity in mapping.values() for shape in entity.get_shapes())) # stable unique via dict shape_tags = [shape.add_to(model.occ) for shape in shapes] # create all shapes before sync @@ -470,9 +471,11 @@ def generate_mesh(model: Model, mapping: Entities) -> None: fragments[shape] = vtags, btags - for name, shape in mapping.items(): - tag = model.addPhysicalGroup(shape.ndims, sorted(shape.select(fragments))) - model.setPhysicalName(dim=shape.ndims, tag=tag, name=name) + for name, entity in mapping.items(): + tag = model.addPhysicalGroup(entity.ndims, sorted(entity.select(fragments))) + model.setPhysicalName(dim=entity.ndims, tag=tag, name=name) + + set_background(model.mesh.field, elemsize, fragments) model.mesh.generate(max(shape.ndims for shape in shapes))