将OpenFOAM中的lduMatrix数据转换为CSC稀疏矩阵格式

将OpenFOAM中的lduMatrix数据转换为CSC稀疏矩阵格式

在OpenFOAM中，lduMatrix是一种用于存储稀疏矩阵的数据结构，主要用于有限体积法离散后的线性系统。要将lduMatrix转换为CSC (Compressed Sparse Column)格式，可以按照以下步骤进行：

理解lduMatrix的存储方式

lduMatrix使用以下主要组件存储矩阵：

• diag(): 对角线元素
• upper(): 上三角非零元素（按行优先顺序存储）
• lower(): 下三角非零元素（按行优先顺序存储）
• lduAddressing: 提供行列索引信息

CSC格式简介

CSC格式由三个数组组成：

1. values: 非零元素值（按列优先顺序）
2. row_indices: 每个非零元素的行索引
3. col_ptrs: 每列起始位置在values和row_indices中的索引

转换步骤

#include "lduMatrix.H"
#include <vector>

void lduToCSC(
    const lduMatrix& matrix,
    std::vector<double>& values,
    std::vector<int>& row_indices,
    std::vector<int>& col_ptrs)
{
    const label nCells = matrix.diag().size();
    const scalarField& diag = matrix.diag();
    const scalarField& upper = matrix.upper();
    const scalarField& lower = matrix.lower();
    const lduAddressing& lduAddr = matrix.lduAddr();
    
    // 首先计算每列的非零元素数
    std::vector<int> nnzPerCol(nCells, 1); // 对角线元素
    
    const labelUList& upperAddr = lduAddr.upperAddr();
    const labelUList& lowerAddr = lduAddr.lowerAddr();
    
    forAll(upperAddr, facei)
    {
        nnzPerCol[upperAddr[facei]]++; // 上三角元素增加对应列计数
        nnzPerCol[lowerAddr[facei]]++; // 下三角元素增加对应列计数
    }
    
    // 构建col_ptrs
    col_ptrs.resize(nCells + 1);
    col_ptrs[0] = 0;
    for (label i = 0; i < nCells; ++i)
    {
        col_ptrs[i+1] = col_ptrs[i] + nnzPerCol[i];
    }
    
    // 分配空间
    values.resize(col_ptrs[nCells]);
    row_indices.resize(col_ptrs[nCells]);
    
    // 临时数组跟踪每列当前插入位置
    std::vector<int> colCurrentPtr(nCells, 0);
    
    // 填充数据
    for (label col = 0; col < nCells; ++col)
    {
        // 首先处理对角线元素
        int pos = col_ptrs[col] + colCurrentPtr[col]++;
        values[pos] = diag[col];
        row_indices[pos] = col;
        
        // 处理上三角元素（列col中的行i > col的元素）
        forAll(upperAddr, facei)
        {
            if (upperAddr[facei] == col)
            {
                int row = lowerAddr[facei];
                int pos = col_ptrs[col] + colCurrentPtr[col]++;
                values[pos] = lower[facei]; // 注意上下三角对应关系
                row_indices[pos] = row;
            }
        }
        
        // 处理下三角元素（列col中的行i < col的元素）
        forAll(lowerAddr, facei)
        {
            if (lowerAddr[facei] == col)
            {
                int row = upperAddr[facei];
                int pos = col_ptrs[col] + colCurrentPtr[col]++;
                values[pos] = upper[facei]; // 注意上下三角对应关系
                row_indices[pos] = row;
            }
        }
    }
}

使用示例

// 假设有一个lduMatrix对象
lduMatrix& matrix = ...;

std::vector<double> values;
std::vector<int> row_indices;
std::vector<int> col_ptrs;

lduToCSC(matrix, values, row_indices, col_ptrs);

// 现在values, row_indices, col_ptrs包含了CSC格式的数据

注意事项

1. OpenFOAM中的upper()和lower()存储顺序需要特别注意，它们对应于面的连接关系。
2. 对称矩阵可以只存储上三角或下三角部分以节省空间。
3. 对于大型矩阵，此转换可能需要优化以提高效率。
4. 如果需要与外部库（如Eigen、SuiteSparse等）交互，确保索引类型匹配（通常是0-based或1-based）。

如果需要更高效的实现或处理特殊情况（如对称矩阵），可能需要对上述代码进行相应调整。