thm_assemble_T_system().cc

    // local ->globe
    cell->get_dof_indices(T_local_dof_indices);
    for (unsigned int i = 0; i < T_dofs_per_cell; ++i) {
      for (unsigned int j = 0; j < T_dofs_per_cell; ++j) {
        auto local_i = T_local_dof_indices[i];
        auto local_j = T_local_dof_indices[j];
        T_mass_matrix.add(local_i, local_j, T_local_mass_matrix(i, j));
        T_stiffness_matrix.add(local_i, local_j,
                               T_local_stiffness_matrix(i, j));
      }
      T_system_rhs(T_local_dof_indices[i]) += T_local_rhs(i);
    }
    T_system_matrix.copy_from(T_mass_matrix);
    T_system_matrix.add(
        theta * time_step,
        T_stiffness_matrix);  // T_mass_matrix +
                              // theta*time_step*T_stiffness_matrix