mesh.NumEntities(codim)

for(const lf::mesh::Entity& entity : mesh.Entities(codim))
	mesh.Index(entity)
   
   entity.SubEntities(subCodim) //get an iterator or convert to base::RandomAccessRange<T>
   entity.NumSubEntities(subCodim)
   entity.Geometry();//A Pointer!

Eigen::SparseMatrix<double> S(N, N)
S.coeffRef(i, j) = value;
S.makeCompressed; //convert to CRS

//Using Mesh
lf::assemble::UniformFEDofHandler dofh(mesh_p, {
    {lf::base::RefEl::kPoint(), 1},
    {lf::base::RefEl::kSegment(), 0},
    {lf::base::RefEl::kTria(), 0},
    {lf::base::RefEl::kQuad(), 0}
  });

Ak = getElementMatrix(entity);
A(dofh(i), dofh(j)) += Ak(i, j);
Phik = getElementVector(entity);
Phi(dofh(i)) += Phi(i);

auto boundary = lf::mesh::utils::flagEntitiesOnBoundary(dofh.Mesh(), 2);
auto boundary_values = [&boundary, &dofh](unsigned int dofId)-> std::tuple<bool, double> {
      return {boundary(dofh.Entity(dofId)), g(.)};
  };
lf::assemble::fix_flagged_solution_components(boundary_values, A, Phi);

//See also docs e.g.
lf::uscalfe::InitEssentialConditionFromFunction(...)

SparseMatrix<double> A(m,n);
std::vector<Triplet> triplets;
triplets.reserve(nnz);

for (...) {
   triplets.push_back(Triplet(i, j, value));

A.setFromTriplets(triplets.begin(), triplets.end());
A.makeCompressed();
   
//OR
A.coeffRef(i,j) += val;
   
//Note! for matrix assembly lf::assemble dont use eigen
lf::COOMatrix<double> A(dofs, dofs);
lf::assemble::Assemble...
Eigen::SparseMatrix<double> A_crs = A.makeSparse();

Eigen::SimplicialLDLT<SparseMatrix> solver;
solver.compute(A); //calculate decomposition for A
if (solver.info() != Eigen::Success)
   throw std::runtime_error("Could not decompose the matrix");
x = solver.solve(b);   //Ax=b

Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>> solver;
solver.analyzePattern(A);
solver.factorize(A);
x = solver.solve(b);

std::tuple<T1, T2...>(val1, val2...);
std::make_tuple(val1, val2,...);
{val1, val2,...} //even quicker, may have be explicit w/ types <...>
std::get<idx>(tuple);