Interpolator¶

def interpolate_all(self):
    """Go through all em forces file and interpolate them"""
    interpolated_results = {}
    for name, values in tqdm(self.src_values.items(), total=len(self.src_values)):
        interpolated, unmapped = self.tree.interpolate(values)
        interpolated_results[name] = {
            "interpolated": interpolated,
            "unmapped": unmapped,
        }
    self.interpolated_results = interpolated_results

def export_to_ansys(self, outdir: Path):
    """Export each interpolated result to an ANSYS format

    Parameters
    ----------
    outdir : Path
        output directory for the interpolated files
    """
    if self.interpolated_results is None:
        raise ValueError(
            "No interpolated results found. Run interpolate_all() first."
        )

    templates = _select_template(self.tree.config.interpolated_load)

    for name, result in self.interpolated_results.items():
        outfile = Path(outdir, f"interpolated_{name}.txt")

        with open(outfile, "w") as f:
            for i, template in enumerate(templates):
                for j, id in enumerate(self.tree.ids):
                    f.write(template.format(id, result["interpolated"][j, i]))

    logging.info(f"ANSYS files exported to {outdir}")

def compute_scalar_integrals(self) -> dict[str, float]:
    """Compute the integrals of the interpolated results over the destination mesh"""
    if self.interpolated_results is None:
        raise ValueError(
            "No interpolated results found. Run interpolate_all() first."
        )

    integrals = {}
    for name, result in self.interpolated_results.items():
        if self.tree.config.num_components > 1:
            raise IncompatibleResultsError(
                "Integral computation only supported for scalar results (num_components=1)."
            )
        if self.volumes is not None:
            integral = np.sum(
                result["interpolated"] * self.volumes[:, np.newaxis], axis=0
            )
        elif self.areas is not None:
            integral = np.sum(
                result["interpolated"] * self.areas[:, np.newaxis], axis=0
            )
        else:
            raise IncompatibleResultsError(
                "No volume or area information found. Cannot compute integrals."
            )
        integrals[name] = integral

    return integrals

def compute_EM_resultants(
    self, pole: np.ndarray | None = None
) -> dict[str, dict[str, float]]:
    """Compute the force and moment resultants of the EM forces and the interpolated
    forces, and their differences

    Parameters
    ----------
    pole : np.ndarray | None, optional
        reference pole for moment calculation, by default None (0,0,0)

    Returns
    -------
    dict[str, dict[str, float]]
        dictionary containing the resultants and their differences for each
        interpolated load
    """
    if self.interpolated_results is None:
        raise ValueError(
            "No interpolated results found. Run interpolate_all() first."
        )
    if pole is None:
        pole = np.array([0.0, 0.0, 0.0])

    resultants = {}
    for name in self.src_values.keys():
        F_EM = self.src_values[name]
        F_Mech = self.interpolated_results[name]["interpolated"]

        R_F_EM = np.sum(F_EM, axis=0)
        R_F_Mech = np.sum(F_Mech, axis=0)

        M_EM = np.cross(self.tree.src_coordinates - pole, F_EM)
        M_Mech = np.cross(self.tree.dest_coordinates - pole, F_Mech)

        R_M_EM = np.sum(M_EM, axis=0)
        R_M_Mech = np.sum(M_Mech, axis=0)

        f_err_comp = np.divide(R_F_EM - R_F_Mech, R_F_EM)
        m_err_comp = np.divide(R_M_EM - R_M_Mech, R_M_EM)

        resultants[name] = {
            "R_F_EM": R_F_EM,
            "R_F_Mech": R_F_Mech,
            "R_M_EM": R_M_EM,
            "R_M_Mech": R_M_Mech,
            "f_err_comp": f_err_comp,
            "m_err_comp": m_err_comp,
            "Unmapped_EM_Force": np.linalg.norm(
                self.interpolated_results[name]["unmapped"]
            ),
        }

    return resultants

def build_vtk_output(self, outdir: Path | None = None):
    """Build the VTK output files for visualization

    Parameters
    ----------
    outdir : Path | None, optional
        output directory for the vtk files, if None, files are not saved, by default None

    """
    if self.interpolated_results is None:
        raise ValueError(
            "No interpolated results found. Run interpolate_all() first."
        )
    # build and dump the vtks
    for name in self.src_values.keys():
        # Destination mesh
        pdata = pv.PolyData(self.tree.dest_coordinates)
        interp_matrix = self.interpolated_results[name]["interpolated"]

        # For scalar data (1 component), flatten to 1D array for better performance
        if interp_matrix.shape[1] == 1:
            pdata["Value"] = interp_matrix.ravel()
        else:
            pdata["Value"] = interp_matrix
            for i in range(interp_matrix.shape[1]):
                pdata[f"Component_{i}"] = interp_matrix[:, i]
        self.dest_vtk[name] = pdata

        # Source mesh
        pdata_em = pv.PolyData(self.tree.src_coordinates)
        interp_matrix_em = self.src_values[name]

        # For scalar data (1 component), flatten to 1D array for better performance
        if interp_matrix_em.shape[1] == 1:
            pdata_em["Value"] = interp_matrix_em.ravel()
        else:
            pdata_em["Value"] = interp_matrix_em
            for i in range(interp_matrix_em.shape[1]):
                pdata_em[f"Component_{i}"] = interp_matrix_em[:, i]

        self.src_vtk[name] = pdata_em

        if outdir is not None:
            outfile = Path(outdir, f"{name}_interpolated.vtk")
            self.dest_vtk[name].save(outfile)
            outfile_em = Path(outdir, f"{name}_src.vtk")
            self.src_vtk[name].save(outfile_em)