Thank you so much for your help! It now goes quite well except a minor glitch, and I guess this mostly likely relates to the former mesh problem. After loading mesh, it showed that the mesh.cmesh.n_coor was 23996, which means there were 23996 nodes in total. Then the domain was created with a warning: bad element orientation, trying to correct... corrected. And then I called omega = domain.create_region('Omega','all'), node_list=omega.entities[:,0]. Ideally the node_list should be in length 23996 and contain all node indices from 0 to 23995. However, len(node_list) = 23989 with several indices missing. Is this because of my bad file file and some nodes were deleted during the domain = FEDomain('domain',mesh) process? The final stiffness matrix K results in shape (23989,23989). I try to utilize this matrix by deleting my original nodes which are not in node_list. Is it appropriate to do so or what should I do to match the stiffness matrix to original nodes?