Running CalculiX solver using WELSIM

CalculiX is a classic open-source, general-purpose finite element solver licensed under the GPL. It supports linear and nonlinear static, dynamic, thermal, and contact analyses, and is natively compatible with the Abaqus INP input file format. From university research to small-to-medium-sized enterprise simulations, CalculiX is widely used. As a general-purpose CAE pre and post-processing software, WELSIM features native pre-processing support for CalculiX: users can either export CalculiX INP input scripts for offline calculation with a single click, or directly call the CalculiX solver within the WELSIM graphical user interface to complete computations, streamlining the workflow of “CAD import -> mesh -> material -> boundary conditions -> solving.” This article outlines the steps for co-solving finite element models utilizing both software applications under the Windows operating system. Software Preparation Download and install the latest version of WELSIM; this blog’s example uses the 2026R3 development version. Press enter or click to view image in full size The CalculiX solver can be downloaded as a pre-compiled package from the official website (www.calculix.de) or GitHub. After extraction, you will find the solver executable files in the bin folder: ccx.exe (serial) and ccx_MT.exe (multi-threaded parallel version). Users do not need to modify any of the extracted CalculiX files. Please avoid placing the CalculiX folder on the system drive (such as the C drive) and refrain from having an excessively long absolute folder path. Path Configuration When executing the CalculiX solver via WELSIM for the first time, you need to briefly configure the CalculiX path. The steps are as follows: Write on Medium 1.Open the WELSIM software, and click Preferences in the menu or toolbar. Press enter or click to view image in full size 2.Locate the Solvers tab in the left-hand list. 3.Find the CalculiX file field on the right, click Browse, and select the extracted ccx_MT.exe or ccx.exe. On Windows, an example path could be: D:\WelSimLLC\executable34\CalculiX-2.23.0-win-x64\bin\ccx_MT.exe. Press enter or click to view image in full size 4.Click OK to save the configuration. 5.During specific finite element analyses, simply set the Solver property of the Study Settings node to CalculiX in the project tree window. Once configuration is complete, you can directly call the CalculiX solver. Conclusion This article describes how to configure the co-solving of finite element models between CalculiX and WELSIM under the Windows operating system. In a Linux environment, the setup and operations are largely identical. The WELSIM installation package does not include the CalculiX solver, so users should download the CalculiX binary program and configure its path themselves as demonstrated in the article. WELSIM currently does not support post-processing for CalculiX results, but reading and displaying CalculiX result files will be supported in future versions. WELSIM and the author are not affiliated with CalculiX or Abaqus, and have no direct relationship with the development teams or organizations of CalculiX and Abaqus. The names of the open-source software referenced here are used solely for technical blog articles and software usage reference.

WELSIM调用有限元求解器CalculiX

CalculiX是一款经典开源通用有限元求解器，基于GPL协议，支持线性/非线性静力、动力学、热、接触等分析，原生兼容Abaqus INP输入文件格式，在高校科研、中小企业仿真中使用广泛。 WELSIM作为一款通用CAE前后处理软件，内置对CalculiX的前处理支持：既可以一键导出CalculiX的INP输入脚本离线计算，也可以在WELSIM图形界面内直接调用ccx求解器完成计算，打通「CAD导入→网格→材料→边界→求解」的流程。 本文介绍Windows操作系统下，两个软件联合求解有限元模型的步骤。 软件准备 下载并安装最新版本的WELSIM，文本使用2026R3开发版作为示例。 图片 CalculiX求解器可以从官网www.calculix.de或者GitHub下载预编译包，解压缩后在bin文件夹下会看到ccx.exe（串行）、ccx_MT.exe（多线程并行版本）的求解器可执行文件。用户无需对CalculiX解压缩后的文件做任何修改。需要注意的是，不要将CalculiX文件夹放置在系统盘如C盘，同时文件夹的绝对路径不宜过长。 配置路径 首次联合WELSIM调用CalculiX求解器时，需要简单配置一下CalculiX的路径。步骤如下： 1. 打开WELSIM软件，点击菜单或者工具栏的Preferences（首选项）。 图片 2. 在左侧列表找到 `Solvers`（求解器）选项卡。 3. 在右侧找到CalculiX file一栏，点击Browse，选中解压好的ccx_MT.exe，或者ccx.exe。Windows下示例路径可以是：D:\WelSimLLC\executable34\CalculiX-2.23.0-win-x64\bin\ccx_MT.exe。 图片 4. 点击OK，保存配置。 5. 具体有限元分析时，在项目结构树窗口中，将Study Settings节点的求解器（Solver）属性设为CalculiX即可。 图片 设置完成，可以直接调用CalculiX求解器。 总结 本文描述了如何在Windows操作系统下配置WELSIM与CalculiX联合求解有限元模型，Linux环境下的设置与操作基本一致。WELSIM安装包不自带CalculiX求解器，需要用户自行下载CalculiX二进制程序并配置路径，这是本文核心内容。 WELSIM还未支持CalculiX结果文件的后处理功能，将会在以后的版本中，支持读取并显示CalculiX的结果文件。 WelSim与作者不隶属于CalculiX, Abaqus。和CalculiX, Abaqus开发团队与机构没有直接关系。这里引用开源软件的名称仅用作技术博客文章与软件使用的参考。

The mesh of engineering simulation 5: STL meshing

STL is a data format for representing outer surface models. It breaks down the surface of a 3D object into a large number of triangular facets. Each triangular facet stores the coordinates of three vertices and a unit normal vector (pointing outward from the model). Curved surfaces are not represented using analytical equations; instead, they are approximated and fitted using a massive number of small facets. Unlike parametric solid CAD data formats, STL only stores the surface and not the internal solid, as it lacks units, colors, materials, textures, CAD feature history, assembly trees, and other information. Press enter or click to view image in full size The data in an STL file must be a closed manifold mesh — otherwise, it cannot be used for subsequent 3D printing or finite element computation. However, the vertices in STL data are stored redundantly, meaning there are no “edges,” no “face ownership,” no sharp angle features, and no concept of inside versus outside. This imposes certain limitations when setting boundary conditions for future finite element analysis, requiring the construction of curved surfaces or sharp edges through parameterization. STL has a vast range of use cases. For example, in 3D printing and additive manufacturing, STL is the de facto standard exchange format. 3D scanners (laser or structured light) output point clouds that are reconstructed into triangular meshes, which can be used for part replication, cultural heritage digitization, and shape inspection. Moreover, it is widely used in geometric inspection, model repair, and visual preview. Press enter or click to view image in full size Mesh Generation Algorithms STL surface mesh data cannot be directly used for finite element computation, as it needs to be remeshed into high-quality volumetric meshes. Typically, there are two classical methods: the Direct Method and the Surface Classification Method. The main difference between these two methods is whether facet clustering is performed on the STL data. The Direct Method has a simpler workflow, performs no sharp-edge recognition, and does not repair geometric defects, but the watertightness of the STL data follows more stringent requirements. The Surface Classfication Method supports more flexible human-computer graphic interaction operations while simultaneously splitting multiple geometric faces to support more complex models. Press enter or click to view image in full size Direct Method The Direct Method involves building a single closed shell directly based on the STL surface mesh. The steps are as follows: 1.Read STL data and build the enclosure topology. Vertex merging (removing duplicate vertices): Each triangle in an STL file stores its vertices independently, resulting in a large number of duplicate coordinate points. Coincident vertices are merged based on a floating-point tolerance to establish a global vertex list. Build edge-triangle adjacency relationships: Iterate through all triangles, extract each edge, and record how many triangular facets share each edge. Topology checking: Identify boundary edges (belonging to only 1 triangle) and internal edges (belonging to 2 triangles); detect non-manifold edges (≥ 3 triangles sharing a single edge). Check if the shell is closed: if isolated boundary edges exist, determine that the shell surface is not closed. Facet normal consistency check: Iterate through the facets and attempt to unify the outward normals. If the normals in the STL file are disorganized, it may lead to errors in ensuing mesh generation. Build shell topology: Assemble the complete boundary triangular facet data structure for subsequent meshing algorithms (such as the Advancing Front method) to read boundaries. Press enter or click to view image in full size 2.Perform 1D edge meshing, 2D surface meshing, and 3D volumetric meshing respectively. For 3D volumetric meshing, the classical Advancing Front method or parallel tetrahedral Delaunay method can be used. Write on Medium 3.Mesh optimization. Optimize the overall mesh by inspecting element quality. Element quality optimization removes flat “thin elements” and improves the convergence of finite element computation. Facet Classification Method In many practical applications, we hope to generate a topology similar to a CAD model to facilitate setting up finite element boundary and other conditions later. In this case, the Facet Classification Method can be used to combine all facets into multiple geometric faces, effectively generating discrete B-Rep data. This method also handles complex or defective STL data better. The steps for the Facet Classification Method are as follows: 1.Read STL data and directly build the shell topology. * Classify triangular facets based on normal angles: If the normal angle between adjacent triangles — namely the dihedral angle — is less than the threshold, they are grouped into the same smooth surface. If it exceeds the threshold, it is determined to be a sharp edge and used to capture corners. * Determine and store the correspondence between edges and triangular facets. * Perform classification and grouping on a large number of discrete STL facets using the surface mesh’s adjacency relationships, dihedral angles, curvature, and non-manifold markers. Triangles belonging to the same smooth geometric surface are categorized into one group, and the dividing boundaries (edges, non-manifold edges) are marked. The output is a set of connected facet regions. Press enter or click to view image in full size 2.Build discrete solid. Based on the categorized facets, reconstruct geometric solids (faces, lines, points) to transform the discrete triangular mesh into a geometric model resembling a B-Rep. The aim is to create a B-Rep geometric entity for each cluster. If high-order elements need to be created, the discrete surface can be fitted into a UV-parameterized surface. This upgrades the pure surface mesh into geometric faces or lines that can be recognized and manipulated by CAD modules. 3.Combine all outer surfaces into a closed surface loop to define a closed shell. Typically, the first surface loop is used as the outer closed shell, and all subsequent surface loops are internal cavities (holes), which are excavated from the outer domain. A 3D solid domain is then built on the basis of this shell to serve future 3D volumetric meshing. 4.Perform 1D edge meshing, 2D surface meshing, and 3D volumetric meshing respectively. (Same as the Direct Method) 5.Mesh optimization. (Same as the Direct Method) Press enter or click to view image in full size Conclusion This article introduced two classical methods for generating finite element meshes from STL models. Unlike CAD models like B-Rep, STL data lacks topological data, considerably complicating mesh generation. In practical engineering, one often encounters defective STL files — such as those that are insufficiently airtight, contain holes, non-manifold edges, self-intersections, or inverted facet normals — which further increase the difficulty of meshing and place high demands on the robustness of the mesher. The algorithms introduced in this text can also be used to repair STL files and reduce various issues that arise during subsequent meshing. Furthermore, these algorithms can also be applied to finite element meshing for other surface models such as OBJ, 3MF, OFF, and PLY.

工程仿真软件中的网格5：STL模型的网格划分

STL是一种外表面模型的数据格式。其原理是把三维物体表面拆成大量三角形面片，每个三角面片存储3个顶点坐标和面法向量（指向模型外侧）。模型曲面不使用解析方程，仅用大量小三角形近似拟合曲面。和参数化的实体 CAD数据格式不同，STL只存表面，不存内部实体；没有单位、颜色、材质、纹理、CAD 特征历史、装配树等信息。 Image STL文件中的数据必须是闭合流形网格，否则无法用于后续的3D打印或有限元计算。同时，STL数据中的顶点重复存储，没有 “边”、没有 “面归属”，没有棱角特征与内部 / 外部的概念，这对后续有限元分析的边界条件设置带来一些限制，需要通过参数化构建曲线表面或棱角边。 STL有着大量的使用场景，如3D 打印 / 增材制造，STL是3D打印事实上的标准交换格式。三维扫描仪（激光 / 结构光）输出点云到重建三角网格，可用于零件复刻、文物数字化、外形检测。此外，在几何检查、模型修复、可视化预览也有广泛应用。 Image 网格划分算法 STL 表面网格数据无法直接用于有限元计算，需要重新划分高质量的体网格。常见经典方法有两种：直接法和面片类聚法。这两种方法的主要区别就是是否对STL数据进行面片类聚，直接法工作流更简单，不做锐边识别，不修复几何缺陷，但对 STL数据水密性要求更严苛；类聚法可以支持更灵活的人机图形交互操作，同时拆分多个几何面后支持更复杂模型。 Image 直接法 直接法即是直接在STL表面网格基础上构建单一封闭壳体。步骤如下： 1. 读取STL数据，构建外壳拓扑结构 顶点合并（去重复顶点）：STL 文件中每个三角独立存储顶点，因此含有大量重复坐标点；按浮点容差合并重合顶点，建立全局顶点列表。 构建边-三角邻接关系。遍历所有三角，提取每条边，记录每条边被几个三角面片共享。 拓扑检查。识别边界边（只属于 1 个三角）；内部边（属于两个三角形）；检出非流形边（≥3 个三角共用一条边）；检查壳体是否闭合：如果存在孤立边界边，判定壳体不封闭。 面片法向一致性校验。遍历面片，尝试统一外法向；如果STL文件的法向混乱，可能会导致后续网格划分错误。 构建外壳拓扑结构。组装完整的边界三角面片数据结构，供后续的网格算法（如前沿推进法）读取边界。 Image 2. 分别进行一维边网格、二维面网格、和三维体网格划分。其中三维网格划分，可以使用经典前沿推进法（Advancing Front），或者并行的四面体Delaunay法。 3. 网格优化。通过检查单元质量，来对整体网格进行优化。单元质量优化，移除扁平 “薄单元”，提升有限元计算收敛性。 面片类聚法 很多实际工作中，我们希望能够生成像CAD模型那样的拓扑结构，以便后续能更加容易设置有限元边界等条件。此时可以使用面片类聚法，将所有面片组合成多个几何面，即生成离散的B-Rep数据，这样一个过程，也能够更好的处理复杂或者有缺陷的STL数据。面片类聚法步骤如下： 1. 读取STL数据，直接构建外壳拓扑结构 按照法向夹角聚类三角面片。相邻三角法向夹角，即二面角（dihedral angle），小于阈值，归为同一个光滑曲面；超过阈值则判定为特征锐边，用于捕捉棱角。 确定、存储边-三角面的对应关系。 基于表面三角网格的邻接关系、二面角、曲率、非流形标记，对大量离散 STL 三角面片做聚类分组；把属于同一个光滑几何面的三角片归为一类，并标记分割边界（锐边、非流形边），输出一组面片连通区域。 Image 2. 建立离散实体图形。 基于分类后的面片，重建几何实体（面、线、点），把离散三角网格变成近似B-Rep的几何模型。方法是为每一个聚类创建一个B-rep 几何实体，如果要创建高阶单元，还可以通过拟合的方法，离散曲面成为UV参数化的曲面；把纯面网格升级成CAD模块可以识别、操作的几何面或线。 3. 将所有外表面组合成闭合表面环，定义一个封闭的壳。通常以第一个面环作为外封闭壳；后续所有传入的面环都是内部空腔（孔洞），从外域中挖去这些空腔。并在这个壳的基础上建立三维实体域。供后续生成 3D 体网格。 4. 分别进行一维边网格、二维面网格、和三维体网格划分。（同直接法） 5. 网格优化。（同直接法） Image 总结 本文介绍了两种划分STL模型有限元网格的经典方法。由于STL数据相比B-Rep等CAD模型少了拓扑信息，给网格划分带来了不少挑战。实际工程中，还会遇到含有缺陷的STL文件，如不够密闭，孔洞、非流形边、自相交、面片法向，更加大了网格划分的难度。对网格划分器的鲁棒性提出较高要求。 文本所介绍算法还可以用于修复STL文件，以减少后续网格划分出现的各种问题。此外，算法也可以用于OBJ，3MF，OFF, PLY等其他表面模型的有限元网格划分。

Mesh of Engineering Simulation 4: Node Renumbering

After finite element mesh generation, a process called node renumbering can be applied. It only changes the order of the nodes’ memory storage, not altering the geometric coordinates. By modifying the node indexing of large sparse finite element matrices, it optimizes non-zero element distribution, bandwidth, profile, and cache hit rates, which consequently reduces matrix bandwidth, memory consumption, and solver runtimes. Node renumbering does not alter the numerical simulation results, but it significantly enhances computational efficiency. Press enter or click to view image in full size From both practical development and application perspectives, this article discusses two established node renumbering algorithms: the Reverse Cuthill-McKee (RCM) algorithm and the Hilbert Space-Filling Curve algorithm. Each method has distinct advantages, making them suitable for different types of matrix computations. Their key features and trade-offs are summarized in the table below: Algorithm | RCM | Hilbert Input Information | Topological connectivity (no coordinates required) | XYZ geometric coordinates (no element connectivity used) Optimization Target | Minimum matrix bandwidth / profile | Spatial locality, CPU/GPU cache Target Solvers | Direct sparse solvers (SPARSELU, MUMPS) | Iterative solvers (CG, GMRES), GPU Multi-Disjoint-Component Performance | Robust | Degrades when components are far apart Computational Overhead | BFS graph traversal (moderate) | Coordinate transformation and large array sorting (typically faster) Matrix Bandwidth | Excellent | Average, often larger Write on Medium Spatial Locality | No guarantee | Excellent Reverse Cuthill-McKee (RCM) Algorithm RCM is a classic graph reordering algorithm. When applied to FEA node renumbering, it reduces the bandwidth and profile of stiffness matrices by pulling non-zero elements closer to the main diagonal. This minimizes memory overhead during solving and boosts the performance of direct sparse solvers. Because RCM relies purely on graph topological connectivity and not spatial coordinates, geographically adjacent nodes are not necessarily indexed sequentially after reordering. Algorithm Flow Graph Construction: Collect participating mesh nodes and construct an adjacency graph (node-node graph). Pseudo-Peripheral Node Search: Locate a pair of nodes with the maximum graph distance to serve as starting points. This minimizes final bandwidth and sets the root for Breadth-First Search (BFS) traversal. Cuthill-McKee BFS Traversal: Perform level-set BFS traversal to obtain an initial ordering sequence. Sequence Reversal: Reverse the sequence to obtain the final RCM node ordering. Implementation Details Graph Filtering & Component Decomposition: Construct the graph using solid mesh nodes only, ignoring isolated nodes. The graph should be undirected, where edges represent vertices sharing at least one element. Decompose the graph into connected components; if a model consists of multiple disconnected segments, process each connected component independently. Pseudo-Peripheral Node Search: The search for pseudo-peripheral nodes utilizes iterative BFS passes rather than a strict computation of the graph diameter. This engineering approximation significantly improves search traversal efficiency. Element Compatibility: RCM performs well on hybrid tetrahedral and hexahedral meshes. It yields optimal results for continuum meshes with high topological connectivity, but performance deteriorates on highly disconnected models with multiple separate components. Press enter or click to view image in full size Hilbert Space-Filling Curve Algorithm Sorting nodes via a 3D Hilbert space-filling curve allows for geographically adjacent nodes to remain close in memory layout. This improves sparse matrix processing by optimizing CPU L1/L2 cache hit rates and memory bandwidth -noticeably improving performance for iterative solvers and large-scale unstructured meshes. The Hilbert method incurs minimal memory overhead; it avoids building massive node-adjacency graphs and requires only node pointers and a small recursion stack. Algorithm Flow Coordinate Extraction: Fetch all mesh nodes and extract their (x, y, z) spatial coordinates. Bounding Box Calculation: Compute the global bounding box for the entire mesh model. Multi-part models can be mapped directly into the same Hilbert space. Octree Recursive Sorting: A 3D Hilbert curve divides a unit cube [0, 1]³ into 2³ = 8 sub-cubes, passing through them sequentially. Node sorting is implemented via an octree recursion. Instead of computing explicit integer Hilbert keys, array partitions are reordered in-place during recursive traversal. Remapping: Generate a permutation array and perform global node and connectivity table remapping (identical to the post-processing phase of RCM). Implementation Details This algorithm avoids generating explicit 64-bit or 128-bit Hilbert integer indices. By applying in-place octree recursive partitioning combined with Gray code transformations directly on node pointer arrays, it prevents large-integer overflow issues. Note that the Hilbert method does not explicitly minimize matrix bandwidth; the resulting bandwidth is often larger than that produced by RCM, generally making its direct solver performance lower compared to RCM. Conclusion Node reordering changes neither the underlying mathematical FEA solution, element/face IDs, nor node coordinates. It optimizes solving efficiency purely by reorganizing memory storage, which maintains identical simulation results. When selecting between these methods: Solver Type: Use RCM for direct solvers; use Hilbert for iterative solvers. Model Topology: For models containing many spatially separate parts, RCM is preferred. If Hilbert must be used on such models, it is recommended to process physical parts in localized sub-domains. Ultimately, implementations of mesh renumbering algorithms may vary, so any approach that can correctly and efficiently reindex nodes while speeding up solver performance is a robust solution.

工程仿真软件中的网格4：节点重新编号

有限元网格划分后，还可以对节点进行重新编号。这种重新编号，不改动几何坐标，只是改变节点数组的存储顺序。从而改变有限元大型稀疏矩阵的节点编号，可以优化稀疏矩阵非零元分布、带宽、轮廓、和缓存命中率，从而降低矩阵带宽，减少内存与求解时间。节点重新编号不会改变仿真计算结果，但是可以提升有限元计算效率。 Image 本文从实际开发与应用的角度，讨论常用的两种节点重编号算法：Reverse Cuthill McKee（RCM）和Hilbert空间填充曲线算法。两种方法各有优点，适用于不同类型的矩阵计算，两种方法的特征与优点如下表所示。 算法名 Reverse Cuthill McKee Hilbert 输入信息 拓扑连通性，不需要坐标 XYZ 几何坐标，不使用单元连通 优化目标 最小矩阵带宽 /profile 空间局部性、CPU/GPU 缓存 适用求解器 直接稀疏求解器 (SPARSELU,MUMPS) 迭代求解器 (CG,GMRES),GPU FE 多分离部件 表现稳健 部件相距远时效果退化 计算开销 BFS 图遍历，中等 坐标变换 + 大数组排序，通常更快 矩阵带宽 优秀 一般，经常更大 空间局部性 无保证 优秀 Reverse Cuthill McKee算法 RCM是一种经典的图重排序算法，用于有限元节点重编号后，可以减小刚度矩阵带宽(bandwidth) / 轮廓 (profile)，即让非零元素尽量靠近主对角线，降低求解内存开销，优化直接稀疏求解器性能。算法只使用纯拓扑的网格连通关系，不使用几何空间坐标，纯图论算法，因此几何邻近的节点在重排以后，编号不一定邻近。 RCM具体算法流程如下： 1. 收集参与重编号的节点，构建节点 顶点邻接图（graph）。 2. 寻找伪外围节点 (pseudo peripheral node)。图上距离最远的一对作为起点，减少最终带宽。作为广度优先搜索（BFS）遍历起点。 3. 执行 Cuthill McKee BFS 分层遍历，得到编号序列。 4. Reverse 反转序列，得到 RCM 编号。 在实际开发RCM算法时，有一些细节，如 • 只对实体网格节点做图，跳过孤立节点。图是无向图，图边代表两个顶点共同属于至少一个单元。对图做连通分量分解，当模型有多个互不相连部件时，每个连通分量独立处理。 • 伪外围节点的寻找使用多次 BFS 迭代，不是严格的全图直径，而是使用工程近似方法，可以提升搜索遍历效率。 • RCM也适合用于四面体、六面体混合网格。对连通性好的连续体网格效果最好。对高度不连通、多分离部件模型，效果退化。 Image Hilbert空间填充曲线算法 使用Hilbert空间填充曲线对节点排序，可以使几何空间上邻近的节点在内存编号上变得接近，从而达到优化稀疏矩阵的目的，尤其是优化 CPU的L1/L2缓存命中率和内存带宽，对迭代求解器和大规模非结构网格收益明显。Hilbert方法内存开销小，不需要构建庞大的顶点邻接图，只存节点的指针，额外内存只有递归栈。 Hilbert具体算法流程如下： 1. 获取模型全部网格节点，提取每个节点的坐标 (x,y,z)。 2. 计算整个网格模型的全局包围盒 （bounding box），对于多部件模型，也可以直接放在同一个 Hilbert空间下。 3. 三维 Hilbert 曲线把单位立方体 [0,1]^3分成 2^3=8 个子立方体，曲线依次穿过这 8 个子立方体。同时，通过八叉树递归算法来实现节点的排序，无需使用额外的Hilbert 整数键值，而是直接在递归过程中用分区重排数组。 4. 生成置换数组，执行全局顶点与连通表重映射，后续流程和RCM 一致。 值得注意的是，这个算法并不生成显式的 64/128 位 Hilbert 整数索引，而是采用原地八叉树递归划分和格雷码变换，原地对顶点指针数组排序，避免大整数溢出问题。Hilbert法不优化矩阵带宽；生成的矩阵带宽经常比 RCMK 更大；直接求解器性能通常不如 RCM。 总结 本文讨论了两种常见的对网格节点重新排序的方法，网格节点重排序不改变有限元数学解，不改变单元 ID、面 ID，也不改变节点坐标，只是通过改变存储顺序提升求解效率，仿真结果完全等同。 这两种方法在选用时，可以根据后续的求解器来选择，直接求解器使用RCM方法，迭代求解器使用Hilbert方法。也可以根据是否多体模型来选择，模型包含大量空间分离独立零件，优先 RCM，如果坚持用 Hilbert，建议按物理实体分片处理。 网格重编号的算法的具体实现并不唯一，只要能够高效且正确地重新编号，并提升求解器效率，就是优异的算法。

The Mesh of Engineering Simulation 3: Element Quality

In simulation software such as Finite Element Analysis (FEA) tools, mesh quality is critical to the computation. High-quality elements improve computational speed, whereas low-quality elements compromise accuracy or may even terminate the computation. Common metrics for evaluating element quality include aspect ratio, Jacobian ratio, maximum corner angle, parallel deviation, skewness, characteristic length, orthogonal quality, and warping factor. This article discusses these common mesh (element) quality parameters from the perspectives of practical application and software development. Aspect Ratio The aspect ratio describes the elongation of an element. Defined as the ratio of the longest edge to the shortest edge, it indicates whether the element is excessively elongated. Elongated elements fail to accurately capture gradients along their short-edge direction, leading to significant errors in bending and boundary layer stress calculations. The ideal aspect ratio is 1, and higher values correspond to poorer element quality. Overall, aspect ratios are useful for checking meshes in boundary layers, thin-walled structures, and stretched regions. Press enter or click to view image in full size For general models, the overall aspect ratio should be less than 20; less than 10 in areas of stress concentration. When the aspect ratio exceeds 20, elements become severely stretched, resulting in poor FEA accuracy and convergence difficulties. In CFD boundary layer regions where the mesh is refined in the normal direction, aspect ratios between 30 and 100 are permissible. The numerical methods for calculating aspect ratios are straightforward. For triangles, the aspect ratio is the ratio of the longest edge to the diameter of the equivalent inscribed circle. For quadrilaterals, it is the ratio of the average lengths of the two pairs of opposite edges: max(L1, L2) / min(L1, L2). Tetrahedra are treated similarly to triangles, with the aspect ratio derived from the ratio of the circumscribed sphere radius (R) to the inscribed sphere radius (r). For hexahedra, the aspect ratio is the ratio of the maximum to the minimum length among the 12 edges. As these methods demonstrate, the aspect ratio is a purely geometric length-based calculation that does not utilize Gaussian points or Jacobian determinants. Jacobian Ratio The Jacobian ratio is the ratio of the minimum Jacobian determinant to the maximum Jacobian determinant within an element. The Jacobian uniformity metric reflects distortion gradients and measures how uniform the mapping is across the element. It is specifically designed to capture nodal distortion and local contraction or expansion in high-order elements, serving as a core quality control indicator for higher-order elements. The value ranges from -1 to 1. A ratio of 1 indicates a uniform, undistorted element, while values ​​closer to zero indicate significant distortion. The empirical standards are as follows: Jacobian Ratio | 0.3–1.0 | 0.1–0.3 | 0–0.1 | ≤ 0 Element Quality | High-quality element; stable solution | Distorted; reduced stress accuracy; prone to oscillation | Severely distorted; not recommended for high-precision calculations | Negative Jacobian / Degenerate; solver error There are two standards for calculating the Jacobian ratio that are based on Gaussian points and corner nodes, respectively. The Gaussian point method aligns with the integration locations of the solver’s stiffness matrix, making it a reliable mesh metric for the finite element method. The corner node method only checks geometric vertices; while computationally efficient, it may fail to detect cases where vertices are positive but integration points are negative, so it is typically used for a precursory geometric screening. Press enter or click to view image in full size For linear elements, the Jacobian determinant is constant throughout the element, so the Jacobian ratio remains 1; this fails to reflect the element’s quality. In this case, the element’s Jacobian determinant should be compared against that of an ideal element in order to yield a Jacobian ratio that is not always 1. Additionally, since the Jacobian determinant is constant within the element, the calculation only needs to be performed once per linear element, without iterating through all nodes or Gaussian points. For nonlinear elements, such as second-order elements with midside nodes or higher-order elements, shape functions are of higher degree. Consequently, the Jacobian matrix varies across different integration points within the element, requiring dynamic calculation by substituting coordinates at each integration point. Maximum Corner Angle The maximum corner angle is an intuitive metric for assessing element quality, capable of identifying elements with extremely obtuse angles or those that are nearly flat (degenerate). The maximum angle within the element is reported, typically in degrees. A larger value indicates more severe element distortion; values ​​closer to the ideal equiangular state are preferable. The optimal angle for a quadrilateral is 90° (a right angle) and a triangle is 60°. Angles exceeding 150° indicate a poor-quality element, and elements with angles exceeding 165° are not recommended for solver. Press enter or click to view image in full size Angles (measured in degrees) can be calculated using the formula: θ = arccos(u·v(|u|·|v|))·180/π, where u and v represent the vectors of the two sides forming the angle. For a triangular element, the angles θ0, θ1, and θ2 at the three vertices are calculated; the maximum value among them — max{θ0, θ1, θ2} — is defined as the maximum corner angle. Similarly, the maximum corner angle of a quadrilateral element is max{θ0, θ1, θ2, θ3}. Write on Medium For 3D elements, the maximum corner angle of the element is determined by identifying the maximum corner angle on each face and then comparing these values. For a tetrahedral element with four triangular faces (each with three interior angles), the maximum of all 12 interior angles is taken. A hexahedral element has six quadrilateral faces (each with four vertex angles), so all 24 angles are evaluated to obtain the maximum corner angle of the hexahedron. Skewness Skewness describes the tilt and distortion of an element, measuring the extent to which a quadrilateral or hexahedron deviates from a perfect rectangle or cube. It is a normalized, comprehensive metric that penalizes both excessively large and excessively small angles. As a key indicator of element quality for finite element solvers, it is one of the most commonly used general metrics. The value ranges from 0 to 1: a value of 0 indicates an ideal, orthogonal, and perfectly regular element, while a value of 1 indicates a completely skewed element unsuitable for finite element analysis. Generally, values ​​below 0.2 indicate high-quality elements; values ​​below 0.5 are acceptable. Values ​​above 0.8 indicate poor element quality, leading to reduced numerical accuracy, convergence difficulties in iterative calculations, and potential stress singularities. Press enter or click to view image in full size The general calculation method for skewness is as follows: S = max{ (θmax − θideal) / (180 − θideal), (θideal − θmin) / θideal } Here, θmax is the maximum corner angle, θmin is the minimum interior angle, and θideal is the ideal interior angle. The ideal internal angle for a triangle is 60° and for a quadrilateral is 90°; the ideal dihedral angle for a tetrahedron is approximately 70.5° and for a hexahedron is 90°. For 3D elements such as tetrahedra and hexahedra, skewness is calculated using the maximum or minimum dihedral angle. This method, also known as the equiangular method, supports all element types, including prisms and pyramids. For triangular and tetrahedral elements, skewness can also be calculated using the equi-volume method. The formula is: S = (V_ideal — V_actual) / V_ideal, where V_ideal is the volume of a standard equilateral element sharing the same circumsphere. For pyramid and wedge elements, equiangular skewness is calculated for all quadrilateral faces and equi-volume skewness for all triangular faces; the maximum value obtained is taken as the skewness. Element Quality Element quality is a comprehensive metric for evaluating mesh quality, with values ​​ranging from 0 to 1. A value of 1 represents a perfect square (2D) or a perfect cube (3D), while a value of 0 indicates zero volume or negative volume (indicating mesh distortion failure). Press enter or click to view image in full size The metric is determined as follows: For 2D quadrilateral or triangular elements, the calculation uses the ratio of the element area to the sum of the squares of the edge lengths: Q = C * (A / ΣL²). The area of ​​a triangular element can be obtained as half the magnitude of the vector cross product. Quadrilateral elements can be divided into two adjacent triangles to calculate their area. For 3D elements, the calculation uses the ratio of the element volume to the square root of the cube of the sum of the squares of the edge lengths: Q = C * (V / (ΣL²)^(3/2)). The volume of a 3D element can be obtained by calculating the determinant of the element’s Jacobian matrix. The value of constant C depends on the element type, as specified below: Element | Triangle | Quadrilateral | Tetrahedron | Hexahedron | Wedge | Pyramid Constant Value | 6.928 | 4.0 | 124.70 | 41.57 | 62.354 | 96 In element quality assessment, a quality value greater than 0.7 indicates an acceptable element. For elements in stress concentration zones, a value greater than 0.8 is desirable. Elements with values ​​below 0.05 are considered severely defective and must be repaired. Conclusion This article discusses the quality of common element types, as calculating and visualizing element quality is an essential feature of modern general-purpose finite element analysis software. These quality metrics are used not only for finite element meshing but also for interactive mesh inspection. By visualizing element quality via contour plots, users can quickly and intuitively assess mesh quality and decide whether to use the current mesh for subsequent calculations. Although this article does not discuss triangular prism elements or wedge elements in detail, their quality metrics are calculated similarly to those of tetrahedral and hexahedral elements and can be determined using the same methods. Other quality metrics — such as parallel deviation, orthogonal quality, and warping factor — are not addressed here.