artorg-unibe-ch / spline_mesher

Meshing tool for homogenised finite elements based on HR-pQCT images
https://artorg-unibe-ch.github.io/spline_mesher/
MIT License
0 stars 0 forks source link

add here calculation of tensor of inertia (partially written in gmsh_mesh_builde... #30

Closed github-actions[bot] closed 1 year ago

github-actions[bot] commented 1 year ago

--> the cortical_ext, cortical_int_sanity will be used to calculate the tensor of inertia

--> return: cort_ext, cort_int with points including ToI

https://github.com/simoneponcioni/spline-mesher/blob/3c9d558608e9e291a52501b64106c622945bd7a1/src/spline_mesher/spline-volume-creation.py#L801


            thickness_tol=180e-3,
            phases=2,
        )

        cortical_v.plot_mhd_slice()
        cortical_ext, cortical_int = cortical_v.volume_splines()
        # Cortical surface sanity check
        cortex = csc.CorticalSanityCheck(MIN_THICKNESS=cortical_v.MIN_THICKNESS,
                                         ext_contour=cortical_ext,
                                         int_contour=cortical_int,
                                         model=cortical_v.filename,
                                         save_plot=False)

        cortical_ext_split = np.array_split(cortical_ext, len(np.unique(cortical_ext[:, 2])))
        cortical_int_split = np.array_split(cortical_int, len(np.unique(cortical_int[:, 2])))
        cortical_int_sanity = np.zeros(np.shape(cortical_int_split))
        for i, _ in enumerate(cortical_ext_split):
            cortical_int_sanity[i][:, :-1] = cortex.cortical_sanity_check(ext_contour=cortical_ext_split[i], int_contour=cortical_int_split[i], iterator=i, show_plots=False)
            cortical_int_sanity[i][:, -1] = cortical_int_split[i][:, -1]
        cortical_int_sanity = cortical_int_sanity.reshape(-1, 3)

        # TODO: add here calculation of tensor of inertia (partially written in gmsh_mesh_builder.py)
        # --> the cortical_ext, cortical_int_sanity will be used to calculate the tensor of inertia
        # --> return: cort_ext, cort_int with points including ToI

        gmsh.initialize()
        gmsh.clear()
        mesher = Mesher(geo_file_path, mesh_file_path)

        cortex_centroid = np.zeros((len(cortical_ext_split), 3))
        cortical_int_sanity_split = np.array_split(cortical_int_sanity, len(np.unique(cortical_int_sanity[:, 2])))
        INTERSECTION_NUMBER = 4
        cortical_ext_centroid = np.zeros((np.shape(cortical_ext_split)[0] + INTERSECTION_NUMBER, np.shape(cortical_ext_split)[1] + INTERSECTION_NUMBER, np.shape(cortical_ext_split)[2]))
        cortical_int_centroid = np.zeros((np.shape(cortical_int_split)[0] + INTERSECTION_NUMBER, np.shape(cortical_int_split)[1] + INTERSECTION_NUMBER, np.shape(cortical_int_split)[2]))
        idx_list_ext = np.zeros((len(cortical_ext_split), 4), dtype=int)
        idx_list_int = np.zeros((len(cortical_ext_split), 4), dtype=int)
        _cnt = 1
        for i, _ in enumerate(cortical_ext_split):
            cortex_centroid[i][:-1] = mesher.polygon_tensor_of_inertia(cortical_ext_split[i], cortical_int_sanity_split[i])
            cortex_centroid[i][-1] = cortical_ext_split[i][0, -1]
            print(f'index:\t{i}')
            cortical_ext_centroid[i], idx_list_ext[i] = mesher.insert_tensor_of_inertia(cortical_ext_split[i], cortex_centroid[i][:-1])
            cortical_int_centroid[i], idx_list_int[i] = mesher.insert_tensor_of_inertia(cortical_int_sanity_split[i], cortex_centroid[i][:-1])

        cortical_ext_msh = np.reshape(cortical_ext_centroid, (-1, 3))
        cortical_int_msh = np.reshape(cortical_int_centroid, (-1, 3))

        cort_ext_pts_tags, cortical_ext_bspline = mesher.gmsh_geometry_formulation(cortical_ext_msh)
        cort_int_pts_tags, cortical_int_bspline = mesher.gmsh_geometry_formulation(cortical_int_msh)

        mesher.factory.synchronize()
        gmsh.fltk.run()
        gmsh.finalize()

if __name__ == "__main__":
    logging.log(logging.WARNING, "Starting meshing script...")
    print("Executing gmsh_spline_mesh.py")
    start = time.time()
    main()
    end = time.time()
    elapsed = end - start
    print("Elapsed time: ", elapsed)
    logging.log(logging.WARNING, "Meshing script finished.")