Spaces:
Sleeping
Sleeping
import Muscat.Containers.ElementsDescription as ED | |
from Muscat.FE.FETools import PrepareFEComputation | |
from Muscat.FE.Fields.FEField import FEField | |
import Muscat.Containers.MeshModificationTools as MMT | |
from Muscat.Containers import MeshGraphTools as MGT | |
from Muscat.Containers.Filters import FilterObjects as FO | |
from Muscat.Containers.Filters import FilterOperators as FOp | |
import Muscat.Containers.MeshInspectionTools as MIT | |
out_fields_names = ['Ux', 'Uy', 'p', 'nut'] | |
in_scalars_names = ['angle_of_attack', 'inlet_velocity'] | |
from Muscat.Containers.NativeTransfer import NativeTransfer | |
import copy | |
import bisect | |
import numpy as np | |
from scipy import spatial | |
import pickle | |
from Muscat.Bridges.CGNSBridge import CGNSToMesh | |
import torch | |
device = torch.device("cpu") | |
num_steps = [2_000] | |
from plaid.containers.dataset import Sample | |
def pretreat_sample(mesh): | |
################### | |
# Compute the skin of the mesh (containing the external boundary and the airfoil boundary) | |
################### | |
MMT.ComputeSkin(mesh, md = 2, inPlace = True) | |
################### | |
# Extract ids of the bar elements corresponding to the airfoil | |
################### | |
ff1 = FO.ElementFilter(zone = lambda p: (-p[:,0]-1.99)) | |
ff2 = FO.ElementFilter(zone = lambda p: (p[:,0]-3.99)) | |
ff3 = FO.ElementFilter(zone = lambda p: (-p[:,1]-1.49)) | |
ff4 = FO.ElementFilter(zone = lambda p: (p[:,1]-1.49)) | |
efAirfoil = FOp.IntersectionFilter(filters=[ff1, ff2, ff3, ff4]) | |
airfoil_ids = efAirfoil.GetIdsToTreat(mesh, "bar2") | |
################### | |
# Preparations | |
################### | |
ext_bound = np.setdiff1d(mesh.elements["bar2"].GetTag("Skin").GetIds(), airfoil_ids) | |
mesh.elements["bar2"].GetTag("External_boundary").SetIds(ext_bound) | |
mesh.elements["bar2"].GetTag("Airfoil").SetIds(airfoil_ids) | |
nfExtBoundary = FO.NodeFilter(eTag = "External_boundary") | |
nodeIndexExtBoundary = nfExtBoundary.GetNodesIndices(mesh) | |
mesh.GetNodalTag("External_boundary").AddToTag(nodeIndexExtBoundary) | |
nfAirfoil = FO.NodeFilter(eTag = "Airfoil") | |
nodeIndexAirfoil = nfAirfoil.GetNodesIndices(mesh) | |
mesh.GetNodalTag("Airfoil").AddToTag(nodeIndexAirfoil) | |
################### | |
# Add node tag for the intrado and extrado | |
################### | |
nfAirfoil = FO.NodeFilter(eTag = "Airfoil") | |
nodeIndexAirfoil = nfAirfoil.GetNodesIndices(mesh) | |
mesh.GetNodalTag("Airfoil").AddToTag(nodeIndexAirfoil) | |
airfoil = ExtractAirfoil(mesh) | |
indices_extrado = airfoil[0][0] | |
indices_intrado = airfoil[0][1] | |
mesh.GetNodalTag("Extrado").AddToTag(indices_extrado) | |
mesh.GetNodalTag("Intrado").AddToTag(indices_intrado) | |
efExtrado = FO.ElementFilter(nTag = "Extrado") | |
efIntrado = FO.ElementFilter(nTag = "Intrado") | |
mesh.elements["bar2"].GetTag("Extrado").SetIds(efExtrado.GetIdsToTreat(mesh, "bar2")) | |
mesh.elements["bar2"].GetTag("Intrado").SetIds(efIntrado.GetIdsToTreat(mesh, "bar2")) | |
def ExtractPathFromMeshOfBars(mesh, startingClosestToPoint, trigo_dir = True): | |
nodeGraph0Airfoild = MGT.ComputeNodeToNodeGraph(mesh, dimensionality=1) | |
nodeGraphAirfoild = [list(nodeGraph0Airfoild[i].keys()) for i in range(nodeGraph0Airfoild.number_of_nodes())] | |
tree = spatial.KDTree(mesh.nodes) | |
_, indicesTrailEdge = tree.query([startingClosestToPoint], k=1) | |
p1init = indicesTrailEdge[0] | |
temp1=mesh.nodes[nodeGraphAirfoild[p1init][0]][1] | |
temp2=mesh.nodes[nodeGraphAirfoild[p1init][1]][1] | |
if trigo_dir: | |
condition = temp1 > temp2 | |
else: | |
condition = temp1 < temp2 | |
if condition: | |
p2 = nodeGraphAirfoild[p1init][0] | |
else: | |
p2 = nodeGraphAirfoild[p1init][1] | |
p1 = p1init | |
path = [p1, p2] | |
while p2 != p1init: | |
p2save = p2 | |
tempArray = np.asarray(nodeGraphAirfoild[p2]) | |
p2 = tempArray[tempArray!=p1][0] | |
p1 = p2save | |
path.append(p2) | |
return path | |
def ExtractAirfoil(mesh): | |
efAirfoil = FO.ElementFilter(elementType=ED.Bar_2, eTag=["Airfoil"]) | |
airfoilMesh = MIT.ExtractElementsByElementFilter(mesh, efAirfoil) | |
path = ExtractPathFromMeshOfBars(airfoilMesh, np.array([1.,0.])) | |
tree = spatial.KDTree(airfoilMesh.nodes[path]) | |
_, indicesLeadEdge = tree.query([[0.,0.]], k=1) | |
indices_extrado = path[:indicesLeadEdge[0]+1] | |
indices_intrado = path[indicesLeadEdge[0]:] | |
indices_airfoil = [indices_extrado, indices_intrado] | |
nodes_extrado = mesh.nodes[indices_extrado] | |
nodes_intrado = mesh.nodes[indices_intrado] | |
nodes_airfoil = [nodes_extrado, nodes_intrado] | |
return indices_airfoil, nodes_airfoil | |
def computeAirfoilCurvAbscissa(airfoil): | |
indices_airfoil = airfoil[0] | |
nodes_airfoil = airfoil[1] | |
curv_abscissa = [] | |
for i in range(2): | |
local_curv_abscissa = np.zeros(len(indices_airfoil[i])) | |
for j in range(1,len(local_curv_abscissa)): | |
local_curv_abscissa[j] = local_curv_abscissa[j-1] + np.linalg.norm(nodes_airfoil[i][j]-nodes_airfoil[i][j-1]) | |
local_curv_abscissa /= local_curv_abscissa[-1] | |
curv_abscissa.append(local_curv_abscissa) | |
return curv_abscissa | |
def MapAirfoil(airfoil_ref, curv_abscissa_ref, curv_abscissa): | |
nodes_airfoil_ref = airfoil_ref[1] | |
dim_nodes = nodes_airfoil_ref[0][0].shape[0] | |
mapped_airfoil = [] | |
for i in range(2): | |
local_mapped_airfoil = np.zeros((len(curv_abscissa[i])-1, dim_nodes)) | |
for j in range(len(curv_abscissa[i])-1): | |
index = max(bisect.bisect_right(curv_abscissa_ref[i], curv_abscissa[i][j]) - 1, 0) | |
a = nodes_airfoil_ref[i][index] | |
b = nodes_airfoil_ref[i][index+1] | |
dl = curv_abscissa[i][j] - curv_abscissa_ref[i][index] | |
dir = (b-a)/np.linalg.norm(b-a) | |
local_mapped_airfoil[j] = a + dl * dir | |
mapped_airfoil.append(local_mapped_airfoil) | |
return mapped_airfoil | |
def GetFieldTransferOpCppStep1(inputField, nbThreads = None): | |
method="Interp/Clamp" | |
nt = NativeTransfer() | |
if nbThreads is not None: | |
nt.SetMaxConcurrency(nbThreads) | |
nt.SetTransferMethod(method) | |
defaultOptions = {"usePointSearch": True, | |
"useElementSearch": False, | |
"useElementSearchFast": False, | |
"useEdgeSearch": True, | |
} | |
options = {} | |
defaultOptions.update(options) | |
dispatch = {"usePointSearch": nt.SetUsePointSearch, | |
"useElementSearch": nt.SetUseElementSearch, | |
"useElementSearchFast": nt.SetUseElementSearchFast, | |
"useEdgeSearch": nt.SetUseEdgeSearch, | |
"DifferentialOperator": nt.SetDifferentialOperator, | |
} | |
for k, v in defaultOptions.items(): | |
if k in dispatch.keys(): | |
dispatch[k](v) | |
else: | |
raise RuntimeError(f"Option {k} not valid") | |
nt.SetSourceFEField(inputField, None) | |
return nt | |
def GetFieldTransferOpCppStep2(nt, targetPoints): | |
nt.SetTargetPoints(targetPoints) | |
nt.Compute() | |
op = nt.GetOperator() | |
status = nt.GetStatus() | |
return op, status | |
def morph_sample(mesh, airfoil_0, curv_abscissa_0): | |
airfoil_1 = ExtractAirfoil(mesh) | |
curv_abscissa_1 = computeAirfoilCurvAbscissa(airfoil_1) | |
mapped_airfoil = MapAirfoil(airfoil_0, curv_abscissa_0, curv_abscissa_1) | |
############################################################## | |
# Compute global target displacement and masks for RBF field morphing | |
############################################################## | |
indices_extrado_to_morph_1 = airfoil_1[0][0][:-1] | |
indices_intrado_to_morph_1 = airfoil_1[0][1][:-1] | |
other_boundary_ids_1 = mesh.GetNodalTag("External_boundary").GetIds() | |
l1 = len(indices_extrado_to_morph_1) | |
l2 = len(indices_intrado_to_morph_1) | |
l3 = len(other_boundary_ids_1) | |
targetDisplacement = np.zeros((l1 + l2 + l3, 2)) | |
targetDisplacementMask = np.zeros((l1 + l2 + l3), dtype = int) | |
targetDisplacement[:l1] = mapped_airfoil[0] - mesh.nodes[indices_extrado_to_morph_1] | |
targetDisplacement[l1:l1+l2] = mapped_airfoil[1] - mesh.nodes[indices_intrado_to_morph_1] | |
targetDisplacement[l1+l2:l1+l2+l3] = np.zeros((l3,2)) | |
targetDisplacementMask[:l1] = indices_extrado_to_morph_1 | |
targetDisplacementMask[l1:l1+l2] = indices_intrado_to_morph_1 | |
targetDisplacementMask[l1+l2:l1+l2+l3] = other_boundary_ids_1 | |
############################################################## | |
# Compute the morphing | |
############################################################## | |
# RBF morphing | |
mesh_nodes = mesh.nodes.copy() | |
morphed_nodes = MMT.Morphing(mesh, targetDisplacement, targetDisplacementMask, radius=None) | |
mesh.nodes = morphed_nodes | |
mesh.nodeFields['X'] = mesh_nodes[:,0] | |
mesh.nodeFields['Y'] = mesh_nodes[:,1] | |
return mesh | |
def project_sample(morphed_mesh, morphed_mesh_0): | |
projected_mesh = copy.deepcopy(morphed_mesh_0) | |
space_, numberings_, _, _ = PrepareFEComputation(morphed_mesh) | |
inputFEField = FEField(name="dummy", mesh=morphed_mesh, space=space_, numbering=numberings_[0]) | |
nt = GetFieldTransferOpCppStep1(inputFEField, 1) | |
FE_interpolation_op, _ = GetFieldTransferOpCppStep2(nt, morphed_mesh_0.nodes) | |
for pfn in out_fields_names + ['X', 'Y']: | |
projected_mesh.nodeFields[pfn] = FE_interpolation_op.dot(morphed_mesh.nodeFields[pfn]) | |
return projected_mesh | |
def pretreat_morph_and_project_mesh(i_sample, dataset, indices, morphed_mesh_0, airfoil_0, curv_abscissa_0): | |
################### | |
## 1) Pretreat data | |
################### | |
sample = Sample.model_validate(pickle.loads(dataset[int(indices[i_sample])]["sample"])) | |
mesh = CGNSToMesh(sample.get_mesh()) | |
pretreat_sample(mesh) | |
################### | |
## 2) Morph data | |
################### | |
morphed_mesh = morph_sample(mesh, airfoil_0, curv_abscissa_0) | |
# print("morphed_mesh =", morphed_mesh) | |
# from Muscat.IO import XdmfWriter as XW | |
# XW.WriteMeshToXdmf("morphed_mesh.xdmf", morphed_mesh) | |
################### | |
## 3) Project data | |
################### | |
projected_mesh = project_sample(morphed_mesh, morphed_mesh_0) | |
# print("projected_mesh =", projected_mesh) | |
# from Muscat.IO import XdmfWriter as XW | |
# XW.WriteMeshToXdmf("projected_mesh.xdmf", projected_mesh) | |
in_scalars = [sample.get_scalar(sn) for sn in in_scalars_names] | |
return [projected_mesh, in_scalars, morphed_mesh.nodes] | |
def infer(dataset, indice, training_data): | |
common_mesh, correlationOperator2c, airfoil_0, curv_abscissa_0, scalar_scalers, scalerX, pca_clouds, pca_fields, y_scalers, kmodels, X_train, y_train = training_data | |
projected_mesh, in_scalars, morphed_mesh_nodes = pretreat_morph_and_project_mesh(indice, dataset = dataset, indices = np.arange(len(dataset)), morphed_mesh_0 = common_mesh, airfoil_0 = airfoil_0, curv_abscissa_0 = curv_abscissa_0) | |
space_, numberings_, _, _ = PrepareFEComputation(common_mesh) | |
inputFEField_0 = FEField(name="dummy", mesh=common_mesh, space=space_, numbering=numberings_[0]) | |
nt_0 = GetFieldTransferOpCppStep1(inputFEField_0, nbThreads = 2) | |
iffo = GetFieldTransferOpCppStep2(nt_0, morphed_mesh_nodes)[0] | |
clouds = np.stack([projected_mesh.nodeFields["X"], projected_mesh.nodeFields["Y"]], axis=1) | |
scalars = np.array(in_scalars) | |
clouds = clouds.reshape(-1, 1) | |
X_pca = np.dot(pca_clouds, correlationOperator2c.dot(clouds)).T | |
X_scalars = scalar_scalers.transform(scalars.reshape(1, -1)) | |
unscaled_X = np.concatenate([X_pca, X_scalars], axis=-1) | |
X = scalerX.transform(unscaled_X) | |
predictions = {} | |
y = [] | |
for i, fn in enumerate(out_fields_names): | |
X_ = torch.tensor(X, dtype=torch.float64).to(device) | |
output_dim = pca_fields[0].shape[0] | |
n_samples = 1 | |
y_pred_i = np.empty((n_samples, output_dim, len(num_steps))) | |
for j in range(output_dim): | |
for k in range(len(num_steps)): | |
y_pred_i[:,j,k] = kmodels[i][j][k](X_).detach().cpu().numpy().squeeze() | |
y.append(y_pred_i) | |
y_pred_i = np.mean(y_pred_i, axis = 2) | |
y_pred_i = y_scalers[i].inverse_transform(y_pred_i) | |
y_pred_common_i = np.dot(y_pred_i, pca_fields[i]).flatten() | |
predictions[fn] = iffo.dot(y_pred_common_i) | |
return predictions | |