Source code for simulai.mesh

# (C) Copyright IBM Corp. 2019, 2020, 2021, 2022.

#    Licensed under the Apache License, Version 2.0 (the "License");
#    you may not use this file except in compliance with the License.
#    You may obtain a copy of the License at

#           http://www.apache.org/licenses/LICENSE-2.0

#     Unless required by applicable law or agreed to in writing, software
#     distributed under the License is distributed on an "AS IS" BASIS,
#     WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#     See the License for the specific language governing permissions and
#     limitations under the License.

import numpy as np
from itertools import product
from collections import OrderedDict
import copy

[docs]class StructuredMesh: def __init__(self, dim_bounds=None, dim_gdl=None, boundary_dim_gdl=None, dim_tags=None): self.n_dim = len(dim_tags) self.dim_tags = dim_tags self.mesh_tags = list() self.elements = OrderedDict() self.boundary_nodes_tags = list() self.boundary_elements = OrderedDict() self.dim_gdl_tags = {tag: gdl for tag, gdl in zip(dim_tags, dim_gdl)} if boundary_dim_gdl: self.boundary_dim_gdl_tags = {tag: gdl for tag, gdl in zip(dim_tags, boundary_dim_gdl)} else: self.boundary_dim_gdl_tags = self.dim_gdl_tags # Constructing axis for bounds, gdl, tag in zip(dim_bounds, dim_gdl, dim_tags): setattr(self, tag, np.linspace(*bounds, gdl + 1)) mesh_matrices = np.meshgrid(*[getattr(self, tag) for tag in dim_tags]) # Constructing the mesh matrices for tag, matrix in zip(dim_tags, mesh_matrices): matrix_tag = tag.capitalize() + '_f' self.mesh_tags.append(matrix_tag) setattr(self, matrix_tag, matrix) # Constructing the mesh elements for tag in self.dim_tags: gdl = self.dim_gdl_tags[tag] domain = getattr(self, tag).copy() subdomains = [domain[i:i+2] for i in range(0, gdl, 1)] setattr(self, tag+'_subdomains', subdomains) for ii, el in enumerate(product(*[getattr(self, tag+'_subdomains') for tag in dim_tags])): self.elements['el_' + str(ii)] = el for bounds, gdl, tag in zip(dim_bounds, boundary_dim_gdl, dim_tags): setattr(self, tag + '_b', np.linspace(*bounds, gdl + 1)) # Constructing the boundaries for dim, dim_tag in enumerate(self.dim_tags): dim_tags = copy.copy(self.dim_tags) lower_bound = getattr(self, dim_tag + '_b').copy()[0] upper_bound = getattr(self, dim_tag + '_b').copy()[-1] lower_boundary = np.meshgrid(np.array([lower_bound]), *self._get_boundaries_curves(but=dim_tag)) upper_boundary = np.meshgrid(np.array([upper_bound]), *self._get_boundaries_curves(but=dim_tag)) dim_tags.remove(dim_tag) setattr(self, dim_tag + '_' + dim_tag + '_b0', lower_boundary[0]) setattr(self, dim_tag + '_' + dim_tag + '_bL', upper_boundary[0]) lower_boundary.pop(0) upper_boundary.pop(0) for ii, otag in enumerate(dim_tags): setattr(self, otag + '_' + dim_tag + '_b0', lower_boundary[ii]) setattr(self, otag + '_' + dim_tag + '_bL', upper_boundary[ii]) self.boundary_nodes_tags.append(otag + '_' + dim_tag + '_b0') self.boundary_nodes_tags.append(otag + '_' + dim_tag + '_bL') # Constructing the boundary elements for ii, bb in enumerate(self.boundary_nodes_tags): boundary_array = getattr(self, bb).copy() gdl = boundary_array.shape[0] tag = bb.split('_')[0] index = self.dim_tags.index(tag) subdomains = {'el' + str(i) + '_' + bb: (boundary_array[i:i + 2], index) for i in range(0, gdl - 1, 1)} self.boundary_elements[bb] = subdomains def _get_boundaries_curves(self, but=None): return [getattr(self, tag + '_b') for tag in self.dim_tags if tag != but]
[docs] def internal_product(self, vector): if isinstance(vector, list): product_list = self.n_dim*(vector,) elif isinstance(vector, tuple): product_list = vector else: raise Exception("The internal product cannot be performed.") return list(product(*product_list))
[docs] def internal_boundary_product(self, vector): vector_ = np.array(vector) if len(vector_.shape) > 1: vector_ = np.array(vector)[:, tag] else: pass if isinstance(vector_, np.ndarray): product_list = (vector_,) elif isinstance(vector_, tuple): product_list = vector_ else: raise Exception("The internal product cannot be performed.") return list(product(*product_list))
[docs] def map_to_element(self, points, reference_interval, el): lower_bound, upper_bound = reference_interval local_el = np.array(el) local_points = np.array(points) dims_max = local_el.max(1) dims_min = local_el.min(1) points_mapped = dims_min + (dims_max - dims_min)*(local_points - lower_bound)\ /(upper_bound - lower_bound) return points_mapped
[docs] def map_to_boundary_element(self, points, reference_interval, el, tag=None): lower_bound, upper_bound = reference_interval local_el = np.array(el) if tag != None : local_points = np.array(points[tag]) else: local_points = np.array(points) dims_max = local_el.max(0) dims_min = local_el.min(0) points_mapped = dims_min + (dims_max - dims_min)*(local_points - lower_bound)\ /(upper_bound - lower_bound) return points_mapped.T