Source code for pyleecan.Methods.Mesh.MeshMat.interface
# -*- coding: utf-8 -*-
from ....Classes.NodeMat import NodeMat
from ....Classes.CellMat import CellMat
from ....Classes.Interpolation import Interpolation
from ....Classes.FPGNSeg import FPGNSeg
from ....Classes.ScalarProductL2 import ScalarProductL2
from ....Classes.RefSegmentP1 import RefSegmentP1
from ....definitions import PACKAGE_NAME
from collections import Counter
import numpy as np
from itertools import combinations
def interface(self, other_mesh):
"""Define a MeshMat object corresponding to the exact intersection between two meshes (nodes must be in both meshes).
Parameters
----------
self : MeshMat
a Mesh object
other_mesh : Mesh
an other Mesh object
Returns
-------
"""
# Dynamic import
new_mesh = self.copy()
new_mesh._is_renum = True
# new_mesh.node = NodeMat()
new_mesh.cell = dict()
for key in self.cell:
# Developer info: IDK if this code works with other than triangle cells. To be checked.
if self.cell[key].nb_node_per_cell == 3: # Triangle case
new_mesh.cell["line"] = CellMat(nb_node_per_cell=2)
interp = Interpolation()
interp.gauss_point = FPGNSeg()
interp.ref_cell = RefSegmentP1()
interp.scalar_product = ScalarProductL2()
new_mesh.cell["line"].interpolation = interp
connect = self.cell[key].get_connectivity()
connect2 = other_mesh.cell[key].get_connectivity()
nodes_tags = np.unique(connect)
other_nodes_tags = np.unique(connect2)
# Find the nodes on the interface (they are in both in and out)
interface_nodes_tags = np.intersect1d(nodes_tags, other_nodes_tags)
nb_interf_nodes = len(interface_nodes_tags)
comb = combinations(range(self.cell[key].nb_node_per_cell), 2)
for duo in list(comb):
col1i = np.mod(duo[0], 3)
col2i = np.mod(duo[1], 3)
col1 = connect[:, col1i]
col2 = connect[:, col2i]
col1_bin = np.zeros(len(col1))
col2_bin = np.zeros(len(col1))
for node in interface_nodes_tags:
Icol1i = np.where(col1 == node)[0]
Icol2i = np.where(col2 == node)[0]
col1_bin[Icol1i] = 1
col2_bin[Icol2i] = 1
# Position in vector where 2 nodes of the same element are on the interface (potential line element)
I_target = np.where(col1_bin + col2_bin == 2)[0]
comb2 = combinations(range(other_mesh.cell[key].nb_node_per_cell), 2)
for duo2 in list(comb2):
col1j = np.mod(duo2[0], 3)
col2j = np.mod(duo2[1], 3)
col12 = connect2[:, col1j]
col22 = connect2[:, col2j]
col12_bin = np.zeros(len(col12))
col22_bin = np.zeros(len(col12))
for node in interface_nodes_tags:
Icol1i = np.where(col12 == node)[0]
Icol2i = np.where(col22 == node)[0]
col12_bin[Icol1i] = 1
col22_bin[Icol2i] = 1
# Same but in the second mesh
I_target2 = np.where(col12_bin + col22_bin == 2)[0]
for itag in I_target2:
e_tag1 = col12[itag]
e_tag2 = col22[itag]
Iline = np.where(
((col1[I_target] == e_tag1) & (col2[I_target] == e_tag2))
| ((col2[I_target] == e_tag1) & (col1[I_target] == e_tag2))
)[0]
if Iline.size != 0:
new_mesh.add_cell([e_tag1, e_tag2], "line")
return new_mesh