09_VisualizeReconstructed.py

import splinepy as sp
import numpy as np
import pandas as pd
import xml.etree.ElementTree as ET
from os import path
import shutil
import subprocess
import os
from splinepy.helpme.integrate import _get_integral_measure, _get_quadrature_information

from miscellaneous.DataPreparation import save_numpy_to_xml
from miscellaneous.error_evaluation import integrate, integrate_multipatch, compute_integral_error, get_solution_vectors, load_geometry, show_multipatch_field

from models.LinearRegression import LinearRegressionModel
from models.RBFInterpolator import RadialBasisRegressionModel
from models.GaussianProcessRegressor import GaussianProcessRegressionModel

VELOCITY_SIZE = 3066
PRESSURE_SIZE = 1399
DEGREE_ELEVATIONS = 1

common_show_options = {
    "cmap": "jet",
    "scalarbar": True,
    "lighting": "off",
    "control_points": False,
    "knots": False,
}

path_to_executable = 'gismo/stokes_recreation'
models = {
    'linregMod': LinearRegressionModel(),
    # 'gaussianProcessMod_DotWhite': GaussianProcessRegressionModel(),
    # 'gaussianProcessMod_RBF': GaussianProcessRegressionModel(),
    # 'gaussianProcessMod_Matern': GaussianProcessRegressionModel(),
    'gaussianProcessMod_RationalQuadratic': GaussianProcessRegressionModel(),
    # 'gaussianProcessMod_ExpSineSquared': GaussianProcessRegressionModel(),
    # 'radialBasisMod_Linear': RadialBasisRegressionModel(),
    'radialBasisMod_thinplatespline': RadialBasisRegressionModel(),
    'radialBasisMod_cubic': RadialBasisRegressionModel(),
    'radialBasisMod_quintic': RadialBasisRegressionModel(),
    # 'radialBasisMod_multiquadric': RadialBasisRegressionModel(),
    # 'radialBasisMod_inversemultiquadric': RadialBasisRegressionModel(),
    # 'radialBasisMod_inversequadratic': RadialBasisRegressionModel(),
    # 'radialBasisMod_gaussian': RadialBasisRegressionModel(),
}

df = pd.read_excel(f'Data/test/parameter_input.xlsx').transpose()
SAMPLES_TEST = df.shape[1]

for model_name in models:
    print(f'Model: \t\t\t{model_name}')

    # Create temp folder
    if not os.path.exists(f'temp_{model_name}'):
        os.mkdir(f'temp_{model_name}')

    # Read POD modes
    data_training = np.load(f'models/trained_models/{model_name}/MetaData.npz')
    # R = data_training['R']
    R = [1, 11, 21]

    for R_ in R:
        print(f'Singular value: \t{R_}')

        model = models[model_name]
        model_velocity, model_pressure = model.load_model(f'{model_name}/velocity_{int(R_)}.pkl', f'{model_name}/pressure_{int(R_)}.pkl')

        # Predict values
        # Predicting test values
        test_velocity_predict, test_pressure_predict = model.predict_test(model_velocity, model_pressure, rescale=True)

        # for sample in range(SAMPLES_TEST):
        # for sample in range(10):
        for sample in range(5):
            save_numpy_to_xml(test_velocity_predict[:,sample], f'temp_{model_name}/predict_velocity_field_{sample}.xml', rows=VELOCITY_SIZE, cols='1')

            save_numpy_to_xml(test_pressure_predict[:,sample], f'temp_{model_name}/predict_pressure_field_{sample}.xml', rows=PRESSURE_SIZE, cols='1')

            # Get all the input files for the evaluation
            geometry_file = f'Data/test/microstructure/index_{sample}.xml'
            velocity_solution = f'Data/test/velocities/velocity_field_{sample}.xml'
            pressure_solution = f'Data/test/pressure/pressure_field_{sample}.xml'
            velocity_predicted = f'temp_{model_name}/predict_velocity_field_{sample}.xml'
            pressure_predicted = f'temp_{model_name}/predict_pressure_field_{sample}.xml'

            print(f'Pressure ranges from {np.min(test_pressure_predict[:,sample])} to {np.max(test_pressure_predict[:,sample])}')
            print(f'Velocity ranges from {np.min(test_velocity_predict[:,sample])} to {np.max(test_velocity_predict[:,sample])}')

            # Create patches for the geometry with gismo script
            program_output = subprocess.Popen(
                [f'./{path_to_executable}', f'-f {geometry_file}', f'-v {velocity_solution}', f'-p {pressure_solution}', f'-w {velocity_predicted}', f'-q {pressure_predicted}', '-e 1', '--no-plot'],
                stdout=subprocess.PIPE,
                text=True
            ).communicate()[0]

            # Recreate the geometry
            PRESSURE_FILE = 'pressure_field_patches.xml'
            VELOCITY_FILE = 'velocity_field_patches.xml'
            PRESSURE_REC_FILE = 'pressure_field_rec_patches.xml'
            VELOCITY_REC_FILE = 'velocity_field_rec_patches.xml'
            GEOMETRY_FILE = geometry_file
            pressure_data = get_solution_vectors(file_path=PRESSURE_FILE)
            velocity_data = get_solution_vectors(file_path=VELOCITY_FILE, two_dimensional=True)
            pressure_rec_data = get_solution_vectors(file_path=PRESSURE_REC_FILE)
            velocity_rec_data = get_solution_vectors(file_path=VELOCITY_REC_FILE, two_dimensional=True)

             # Load the geometry
            microstructure, ms_vel = load_geometry(GEOMETRY_FILE, degree_elevations=DEGREE_ELEVATIONS)
            
            # Show pressure and velocity field
            mp_pressure = show_multipatch_field(microstructure, pressure_data, data_name="pressure")
            mp_pressure.show()
            # sp.io.default.export('test.svg', mp_pressure)
            mp_velocity = show_multipatch_field(ms_vel, velocity_data, data_name="velocity")
            mp_velocity.show()

    shutil.rmtree(f'temp_{model_name}')
    os.remove(PRESSURE_FILE)
    os.remove(VELOCITY_FILE)
    os.remove(PRESSURE_REC_FILE)
    os.remove(VELOCITY_REC_FILE)