I am using the Eigen library to solve a FEM problem in which I use several similar sparse matrices for the different kinds of derivatives I calculate. To build the sparse matrix to solve the system I would like to use the comma initializer, but that is not supported for sparse matrices (https://eigen.tuxfamily.org/dox/group__TutorialSparse.html and https://stackoverflow.com/a/43532618/15811117). In that answer, Henri Menke does suggest the Eigen unsupported Kronecker product, though I don't think that will work here.
With an NxM mesh, I calculate several derivative matrices (which are sparse), Dx2, Dxz, and Dz2 of dimension N*M x N*M. My solver matrix is as follows, (K's are just constants, and 0 is just filler--a zero matrix, N*M x N*M)
K3*Dx2 + K1*Dz2 K2*Dxz 0 0 0
K2*Dxz K3*Dz2 + K1*Dx2 0 0 0
0 0 K1*(Dx2+Dz2) 0 0
0 0 0 K3*Dx2 + K1*Dz2 K2*Dxz
0 0 0 K2*Dxz K3*Dz2 + K1*Dx2
As of right now, I build it by converting the sparse into dense, using the comma initializer, and then converting back to sparse by .sparseView(). The only problem with this is that if the mesh is much larger than 16x16, I get a std::bad_alloc error (which is reasonably expected, since it creates a huge matrix: 5*N*M x 5*N*M), and I need a mesh at least 100x100 (and there will not always be as many 0 matrices as currently shown, though it will remain a sparse matrix).
What is the best way to build this matrix? I have thought of trying to use Triplets to do so, similar to this (Convert an Eigen matrix to Triplet form C++).
Edit: I don't think that using Triplets will be that easy (though I am currently attempting it) because I edit several of the values in each D_ _ matrix for boundary conditions...unless there is a way to find a Triplet in an std::vector that has given row and column number?
Update: Both the Triplet and the Kronecker product methods are feasible, but I have run into several problems:
With the Kronecker product, I get the following error (about a hundred times or so) at compilation:
.../eigen/unsupported/Eigen/src/KroneckerProduct/KroneckerTensorProduct.h:225:97: error: no type named ‘ReturnType’ in ‘struct Eigen::ScalarBinaryOpTraits<int, std::complex<double>, Eigen::internal::scalar_product_op<int, std::complex<double> > >’
With the Triplet method, I defined some operators and a function to 'shift' the matrix to the right place (basically what the Kronecker product is doing). I just have to figure out an index out of bounds issue.
The problem comes from the fact that the matrices I used in the Kronecker product were of type SparseMatrix<int> and SparseMatrix<complex<double>>. There are no operator overloads for int and complex<double>, so by making both matrices of compatible types, the code compiles and runs as expected.