如何有效地组合有限元稀疏矩阵

How to efficiently assemble a FEM sparse matrix

本文关键字：有限元 有效地 组合      更新时间：2023-10-16

全部处理，感谢您抽出时间阅读我的问题。我使用的是Eigeng3.3.4(http://eigen.tuxfamily.org/index.php?title=Main_Page)编写一些FEM代码。

我阅读了Eigen3.3.4文档，在这个网站上(http://eigen.tuxfamily.org/dox/TopicFunctionTakingEigenTypes.html)，上面写着我们应该使用CCD_ 1来避免额外拷贝并获得高性能。

所以在我的FEM代码中，对于稀疏矩阵组装部分，假设函数为：

FormFE(const Ref<VectorXd> &U,const Ref<VectorXd> &V,
Ref<SparseMatrix<double> > AMATRIX,Ref<VectorXd> RHS)

其中U表示位移，V表示速度项。AMIX是我的稀疏矩阵，RHS是残差项。

然后，我尝试在组装之前首先初始化我的AMIX(我有一个包含所有非零元素及其值的tripleList(我将值设置为零进行初始化))所以我尝试了：

AMATRIX.setFromTriplets(ZeroTripList.begin(),ZeroTripList.end());

但我有一个错误：

class Eigen::Ref<Eigen::SparseMatrix<double, 0, int> >’ has no member named ‘setFromTriplets

那么我该如何解决这个问题呢？

我的解决方案之一是使用：

FormFE(const Ref<VectorXd> &U,const Ref<VectorXd> &V,
SparseMatrix<double> &AMATRIX,Ref<VectorXd> RHS)

这很有效，但我不确定它是否有效。我不太擅长cpp:P

事实上，我的问题是：

  1. 如何有效地使用Eigen(尤其是用于FEM计算)，我在每个FEM相关函数中几乎处处使用Eigen的VectorXd和MatrixXd
  2. 如何有效地组装稀疏矩阵
  3. 有可能为FEM组装进行一些OpenMP并行化吗
  4. 欢迎对基于C++的FEM编码提出任何有用的建议(库推荐或任何有用的想法)

谢谢。顺致敬意，

是的，在这里传递SparseMatrix<double> &是正确的做法。Ref<SparseMatrix>的目的是传递集合到SparseMatrix的集合对象，如子稀疏矩阵、Map<SparseMatrix>。。。

使用setFromTriplets也是确保获得良好性能的正确做法。如果操作得当(即，正确调用reserve和正确的插入顺序)，使用mat.insert(i,j) = val;直接插入元素可能会快x2。但如果你弄错了，它也可能慢x100倍。。。请参阅文档。使用SparseMatrix::insert，也可以使用OpenMP填充矩阵，但这需要更加小心和严格，这里有一个典型的模式：

int n_cols = ??, n_rows = ??;
std::vector<int> nnz_per_col(n_cols);
// set each nnz_per_col[j] to the exact number
// of non-zero entries in the j-th column (or more, but NOT less)
SparseMatrix<double> mat(n_rows, n_cols);
#pragma omp parallel for
for(int j=0; j<cols; ++j) {
for each non zero entry i in the j-th column {
// preferably with increasing i
double val_i_j = ...;
mat.insert(i,j) = val_i_j;
}
}

当然，如果这对你来说更容易的话，你也可以按行工作。在这种情况下，使用SparseMatrix<double,RowMajor>。当然，您可以调整这个模式来处理列/行等块。

如果对于汇编，你需要使用一些密集的矩阵/向量，那么我相信它们非常小，大小固定。然后，与其使用MatrixXd/VectorXd，不如使用静态分配的Ref<MatrixBase>0和Matrix<double,N,1>类型。这将防止大量内存分配/释放。

最后，最重要的建议：如果您关心性能，请不要忘记评测您的代码，然后再调查优化代码的时间和精力。此外，始终打开编译器优化的测试台/配置文件。

相关文章：