diff --git a/src/bgem/gmsh/field.py b/src/bgem/gmsh/field.py index 9190847..95b7218 100644 --- a/src/bgem/gmsh/field.py +++ b/src/bgem/gmsh/field.py @@ -347,6 +347,21 @@ def construct(self, model): # creata MathEval +# def constant(value): +# """ +# Define constant field. +# Unfortunately much slower then the Box field alternative. +# Seems that MathExpr evaluation time is proportional to the string length, i.e. it is parsed every time. +# :param value: +# :return: +# """ +# return FieldExpr(str(value), []) + +def constant(value): + """ + Make a field with constant value = value. + """ + return box((-1, -1, -1), (1, 1, 1), v_in=value, v_out=value) """ Predefined coordinate fields. @@ -414,11 +429,6 @@ def box(pt_min: Point, pt_max: Point, v_in:float, v_out:float=1e300) -> Field: -def constant(value): - """ - Make a field with constant value = value. - """ - return box((-1, -1, -1), (1, 1, 1), v_in=value, v_out=value) @@ -549,30 +559,41 @@ def minimum(*fields: Field) -> Field: # Structured -def geometric(field:Field, a:Tuple[float, float], b:Tuple[float, float]) -> Field: + + +def linear(x_field:Field, a:Tuple[float, float], b:Tuple[float, float]) -> Field: """ - Geometric interpolation between points a = (x0, y0) and b = (x1, y1). + Linear transformation of the scalar field given by two points. + Prescribed mapping a[0] to a[1], and b[0] to b[1]. """ - x = field x0, y0 = a x1, y1 = b - l0, l1 = math.log(y0), math.log(y1) - a = (l1 - l0) / (x1 - x0) - x_avg = (l1 * x0 - l0 * x1) / (l1 - l0) - return current_module.exp(a * (x - x_avg)) + a = (y1 - y0) / (x1 - x0) + #x_avg = (y1 * x0 - y0 * x1) / (y1 - y0) + x_avg = (x0 + x1) / 2.0 + y_avg = (y0 + y1) / 2.0 + + if a > 1e10: + # Following formula could provide more stable results + # but needs one more scalar addition + # that makes the test 3D meshing about 2% slower. + return a * (x_field - x_avg) + y_avg + else: + # Just two field operations. + add = - (a * x0 - y0) + return a * x_field + add -def linear(field:Field, a:Tuple[float, float], b:Tuple[float, float]) -> Field: +def geometric(field:Field, a:Tuple[float, float], b:Tuple[float, float]) -> Field: """ - Linear transformation of the scalar field given by two points. - Prescribed mapping a[0] to a[1], and b[0] to b[1]. + Geometric interpolation between points a = (x0, y0) and b = (x1, y1). + """ x = field x0, y0 = a x1, y1 = b - a = (y1 - y0) / (x1 - x0) - x_avg = (y1 * x0 - y0 * x1) / (y1 - y0) - return a * (x - x_avg) + l0, l1 = math.log(y0), math.log(y1) + return current_module.exp(linear(field, (x0, l0), (x1, l1))) def polynomial(field:Field, a:Tuple[float, float], b:Tuple[float, float], q: float = 1) -> Field: """ @@ -583,9 +604,7 @@ def polynomial(field:Field, a:Tuple[float, float], b:Tuple[float, float], q: flo x0, y0 = a x1, y1 = b l0, l1 = y0 ** (1/q), y1 ** (1/q) - a = (l1 - l0) / (x1 - x0) - x_avg = (l1 * x0 - l0 * x1) / (l1 - l0) - return (a * (x - x_avg)) ** q + return linear(field, (x0, l0), (x1, l1)) ** q def threshold(field:Field, lower_bound:Tuple[float, float], upper_bound:Tuple[float, float]=None, sigmoid:bool=False) -> Field: diff --git a/src/bgem/gmsh/gmsh.py b/src/bgem/gmsh/gmsh.py index 412891c..ce2e73b 100644 --- a/src/bgem/gmsh/gmsh.py +++ b/src/bgem/gmsh/gmsh.py @@ -7,6 +7,8 @@ import gmsh import re import warnings +import inspect +import pathlib from bgem.gmsh import gmsh_exceptions from bgem.gmsh import options as gmsh_options @@ -485,6 +487,9 @@ def import_shapes(self, fileName, highestDimOnly=True): :param highestDimOnly: """ + file_path = pathlib.Path(fileName) + if not file_path.exists(): + raise FileNotFoundError(f"Error: The file '{file_path}' does not exist!") shapes = self.model.importShapes(fileName, highestDimOnly=highestDimOnly) self._need_synchronize = True return ObjectSet(self, shapes, [Region.default_region[dim] for dim, _ in shapes]) @@ -659,6 +664,80 @@ def make_mesh(self, objects: List['ObjectSet'], dim=3, eliminate=True) -> None: if bad_entities: print("Bad entities:", bad_entities) + def get_elements(self, obj_set: 'ObjectSet'): + """ + Get a (dim, tag) dict -> `elementTypes', `elementTags', `nodeTags' + :param obj_set: + :return: + """ + return { + (dim, tag): gmsh.model.mesh.getElements(dim, tag) + for dim, tag in obj_set.dimtags + } + + def mesh_new_region_from_elements(self, dim, region_name, el_tags, el_coords): + el_type = [15, 1, 2, 4][dim] + i_entity = self.model.addDiscreteEntity(dim, -1, []) + self.model.mesh.addElementsByType(i_entity, el_type, el_tags, el_coords) + + # def mesh_level_field_regions(self, obj_set:'ObjectSet', level_fn=None, bins:List[float] = None, region: Union[str, List[str]] = None): + # """ + # Assign region to the elements of the given model entities (dimtags). + # level_fn is given, the regions in form "_" are assigned to the elements + # on which the level_fn is in the bin [x_i, x_{i+1}]. The level_fn is evaluated in nodes, + # elements on the boundary are assigned to the bin for which the nodes have larger total weight. + # Weight of node X is distance of level_fn(X) from the bin borders. + # :param model_entities: + # :param level_fn: a scalar level-set function on set of 3D points, Nx3 -> N + # :param bins: [x_0, x_1, .... ,x_n] define n bins + # default single bin [0,1] + # :param region: str + # default: region = "levelset" + # for List of regions, the ragions assigned to individual bins + # :return: + # + # TODO: It is more tricky to implement as GMSH could assign a single PhysicalGroup to an entity. + # See tutorial X2 from GMSH manual. + # So one has to create a new discrete entity for each discrete level bin. + # gmsh.model.addDiscreteEntity(dim, tag) + # Boundary entities may be specified to allo GMSH perform correct topology operations. + # No need in our case as we will just write down the mesh. + # Then we could possibly assign elements: + # gmsh.model.mesh.addElementsByType(tag, elementType, elementTags, nodeTags): + # """ + # entity_dict = self.get_elements(obj_set) + # node_dict = {} + # for dim,tag in entity_dict.keys(): + # node_ids, node_coords = self.model.getNodes(dim, tag) + # node_dict.update(dict(zip(node_ids, node_coords))) + # unique_ids, unique_coords = zip((k, v) for k, v in node_dict.items()) + # unique_levels = level_fn(unique_coords) + # node_levels = dict(unique_ids, unique_levels) + # + # self.model.model.addPhysicalGroup + # + # + # def mean_level(nodes): + # return np.mean(node_levels[node] for node in nodes) + # + # def bin(level): + # return 1 + # def reg_el_tags(el_tags, el_nodes): + # el_regions = [bin(mean_level(nodes)) for nodes in el_nodes] + # zip(el_tags, el_regions) + # + # + # + # el_data = zip(el_tuple for el_tuple in entity_dict.values()) + # concatenate = lambda x : list(itertools.chain.from_iterable(x)) + # el_ids, el_tags, el_nodes = [concatenate(x) for x in el_data] + # self.mesh_new_region_from_elements(dim, region, el_tags, el_nodes) + # + # + # unique_bins = bin(unique_levels) + # {n_id, for n_id, coord in node_dict.items() + + def write_brep(self, filename=None): self.synchronize() if filename is None: @@ -743,7 +822,6 @@ def group(self, *obj_list: Union['ObjectSet', List['ObjectSet']]) -> 'ObjectSet' - class ObjectSet: default_mesh_step = 0 @@ -793,6 +871,18 @@ def group(*obj_list: 'ObjectSet') -> 'ObjectSet': g.mesh_step_size = mesh_step_size return g + @property + def attr(self): + """ + A dynamic dataclass of attribute maps. + - single attribute : defined name and aggregation operation, own global dimtag map + - global map assigning attribute dicts to dimtags + or track dimtags through operations with associated attributes. + However region and mesh_Step tracking didn't quite work. + - + :return: + """ + pass def __init__(self, factory: 'GeometryOCC', dim_tags: List[DimTag], regions: List[Region]) -> None: self.factory = factory @@ -818,6 +908,10 @@ def tags(self): def size(self): return len(self.dim_tags) + def min_dim(self): + dims = [d for d, t in self.dim_tags] + return min(dims, default=1) + def max_dim(self): dims = [d for d, t in self.dim_tags] return max(dims, default=1) @@ -850,6 +944,59 @@ def modify_regions(self, format: str): self.regions = regions return self + def dt_equal(self, other: 'ObjectSet'): + """ Tests two ObjectSets equality over dimtags. + Dimtags does not need to be sorted. + """ + assert other + return sorted(self.dim_tags) == sorted(other.dim_tags) + + def dt_intersection(self, *obj_list: 'ObjectSet') -> 'ObjectSet': + """ + Create intersection self and given list of ObjectSets (its dimtags) over dimtags. + :param obj_list: List of ObjectSets to be intersected. + :return: new ObjectSet + """ + assert obj_list + result = ObjectSet(factory=self.factory, dim_tags=[], regions=[]) + for item in obj_list: + if isinstance(item, ObjectSet): + for dt in item.dim_tags: + try: + idx = self.dim_tags.index(dt) + result.dim_tags.append(self.dim_tags[idx]) + result.regions.append(self.regions[idx]) + result.mesh_step_size.append(self.mesh_step_size[idx]) + except ValueError as err: + message = "Intersection skip dimtag: {} - not in self object.".format(dt) + # print(message) + else: + raise Exception(f"group: Wrong argument of type {type(item)}, expecting ObjectSet..") + return result + + def dt_drop(self, *obj_list: 'ObjectSet') -> 'ObjectSet': + """ + Drop any number of ObjectSets (its dimtags) from the self object. + If dimtags not found, it is skipped. + :param obj_list: List of ObjectSets which dimtags should be dropped. + :return: self + """ + assert obj_list + for item in obj_list: + if isinstance(item, ObjectSet): + for dt in item.dim_tags: + try: + idx = self.dim_tags.index(dt) + self.dim_tags.pop(idx) + self.regions.pop(idx) + self.mesh_step_size.pop(idx) + except ValueError as err: + message = "Drop skip dimtag: {} - not in self object.".format(dt) + print(message) + else: + raise Exception(f"group: Wrong argument of type {type(item)}, expecting ObjectSet..") + return self + def translate(self, vector): self.factory.model.translate(self.dim_tags, *vector) self.factory._need_synchronize = True @@ -909,16 +1056,116 @@ def revolve(self, center, axis, angle, numElements=[], heights=[], recombine=Fal # split the Objectset by dimtags return all_obj.split_by_dimension() + def fillet(self, curve_obj, radii, removeVolume=True): + """ + gmsh.model.occ.fillet(volumeTags, curveTags, radii, removeVolume=True) + + Fillet the volumes `volumeTags' on the curves `curveTags' with `radii'. + It has effect on the edge curves of the volume rounding them with given radius. + + The `radii' vector can either contain a single radius, as many + radii as `curveTags', or twice as many as `curveTags' (in which case + different radii are provided for the begin and end points of the curves). + Return the filleted entities in `outDimTags' as a vector of (dim, tag) + pairs. Remove the original volume if `removeVolume' is set. + + Return `outDimTags'. + + Types: + - `volumeTags': vector of integers + - `curveTags': vector of integers + - `radii': vector of doubles + - `outDimTags': vector of pairs of integers + - `removeVolume': boolean + """ + assert self.min_dim() == self.max_dim() == 3 + volume_tags = self.tags + assert curve_obj.min_dim() == curve_obj.max_dim() == 1 + curve_tags = curve_obj.tags + radii = np.atleast_1d(radii) + try: + outDimTags = self.factory.model.fillet(volume_tags, curve_tags, radii, removeVolume) + except ValueError as err: + message = "\nFillet failed!\ndimtags: {}".format(str(self.dim_tags[:10])) + gerr = gmsh_exceptions.BoolOperationError(message) + self._raise_gmsh_exception(gerr, err) + + regions = [Region.default_region[dim] for dim, tag in outDimTags] + all_obj = ObjectSet(self.factory, outDimTags, regions) + + self.factory._need_synchronize = True + # split the Objectset by dimtags + return all_obj + + def chamfer(volumeTags, curveTags, surfaceTags, distances, removeVolume=True): + """ + gmsh.model.occ.chamfer(volumeTags, curveTags, surfaceTags, distances, removeVolume=True) + + Chamfer the volumes `volumeTags' on the curves `curveTags' with distances + `distances' measured on surfaces `surfaceTags'. The `distances' vector can + either contain a single distance, as many distances as `curveTags' and + `surfaceTags', or twice as many as `curveTags' and `surfaceTags' (in which + case the first in each pair is measured on the corresponding surface in + `surfaceTags', the other on the other adjacent surface). Return the + chamfered entities in `outDimTags'. Remove the original volume if + `removeVolume' is set. + + Return `outDimTags'. + + Types: + - `volumeTags': vector of integers + - `curveTags': vector of integers + - `surfaceTags': vector of integers + - `distances': vector of doubles + - `outDimTags': vector of pairs of integers + - `removeVolume': boolean + """ + pass + + def get_mass(self): + """ + Total mass of the object dim tags, no sense for mixed dimensions + :return: + """ + return sum(self.factory.model.getMass(dim, tag) for dim, tag in self.dim_tags) + def copy(self) -> 'ObjectSet': + dim_lists = self.split_by_dimension() + copy_dimensions = [ + dim_obj._copy_simple() + for dim_obj in dim_lists + ] + return ObjectSet.group(*copy_dimensions) + def _copy_simple(self) -> 'ObjectSet': """ Problem: gmsh.model.occ.copy fails to copy boundary dimtags. + Problem: it seems to copy only highest dimension tags even if the lower dim tags are not part of the boundary. """ copy_tags = self.factory.model.copy(self.dim_tags) self.factory._need_synchronize = True - copy_obj = ObjectSet(self.factory, copy_tags, self.regions) + copy_regions = self.regions + copy_obj = ObjectSet(self.factory, copy_tags, copy_regions) copy_obj.mesh_step_size = self.mesh_step_size.copy() return copy_obj + def unique(self): + """ + Unique dim_tag list. + :return: + """ + seen = set() + dimtags = [] + regions = [] + + # remove duplicate dimtags + for dt, r in zip(self.dim_tags, self.regions): + if dt not in seen: + seen.add(dt) + dimtags.append(dt) + regions.append(r) + + return ObjectSet(self.factory, dimtags, regions) + def get_boundary(self, combined=False): """ Get the boundary of the model entities dimTags. @@ -1228,6 +1475,40 @@ def remove_small_mass(self, mass_limit): self.regions = regions return self +def _transfrom_op(op, obj_list, *args, **kwargs): + if isinstance(obj_list, (list, tuple, set)): + type_constructor = type(obj_list) + return type_constructor( + _transfrom_op(op, item, *args, **kwargs) + for item in obj_list) + elif isinstance(obj_list, dict): + return { + k: _transfrom_op(op, item, *args, **kwargs) + for k, item in obj_list + } + else: + return op(obj_list, *args, **kwargs) + +def _apply_recursive(fn): + def recursive_fn(data, *args, **kwargs): + return _transfrom_op(fn, data, *args, **kwargs) + return recursive_fn +# Apply transformop to all ObjectSet methods + +for name, method in inspect.getmembers(ObjectSet, predicate=inspect.isfunction): + globals()[name] = _apply_recursive(method) + + +# def translate(obj_list, vector): +# return _transfrom_op(ObjectSet.translate, obj_list, vector) +# +# def rotate(obj_list, axis, angle, center=[0, 0, 0]): +# return _transfrom_op(ObjectSet.rotate, obj_list, axis, angle, center) +# +# def scale(obj_list, scale_vector, center=[0, 0, 0]): +# return _transfrom_op(ObjectSet.scale, obj_list, scale_vector, center) + + def rotation_matrix(axis, theta): """