diff --git a/include/eqlib/Log.h b/include/eqlib/Log.h
index bf190a69..0c26ad0d 100644
--- a/include/eqlib/Log.h
+++ b/include/eqlib/Log.h
@@ -63,14 +63,14 @@ class Log {
     template <class... TArgs>
     static void task_info(std::string text, TArgs&&... args)
     {
-        std::string message = "\033[38;5;33m" + text + "\033[0m";
+        std::string message = "\033[37m" + text + "\033[0m";
         info(2, message.c_str(), std::forward<TArgs>(args)...);
     }
 
     template <class... TArgs>
     static void task_step(std::string text, TArgs&&... args)
     {
-        std::string message = "\033[38;5;8m" + text + "\033[0m";
+        std::string message = "\033[33m" + text + "\033[0m";
         info(3, message.c_str(), std::forward<TArgs>(args)...);
     }
 
@@ -109,7 +109,7 @@ class Log {
     template <class... TArgs>
     static void warn(std::string text, TArgs&&... args)
     {
-        std::string message = "\033[38;5;208m" + text + "\033[0m";
+        std::string message = "\033[35m" + text + "\033[0m";
         s_console->warn(message.c_str(), std::forward<TArgs>(args)...);
     }
 
diff --git a/include/eqlib/Problem.h b/include/eqlib/Problem.h
index 9dd07507..cb854890 100644
--- a/include/eqlib/Problem.h
+++ b/include/eqlib/Problem.h
@@ -78,6 +78,8 @@ class Problem {
     std::vector<std::vector<Index>> m_element_g_equation_indices;
     std::vector<std::vector<Index>> m_element_g_variable_indices;
 
+    std::vector<std::vector<index>> m_element_f_variable_indices_hi;
+
     SparseStructure<double, int, true> m_structure_dg;
     SparseStructure<double, int, true> m_structure_hm;
 
@@ -407,12 +409,34 @@ class Problem {
 
         Log::task_info("The problem occupies {} MB", m_data.values().size() * 8.0 / 1'024 / 1'024);
 
+        Log::task_step("Initialize linear solver...");
+
         #ifdef EQLIB_USE_MKL
         m_linear_solver = new_<PardisoLDLT>();
         #else
         m_linear_solver = new_<SimplicialLDLT>();
         #endif
 
+        Log::task_step("Initialize element boundaries...");
+
+        m_element_f_variable_indices_hi.resize(nb_elements_f);
+
+        for (index i = 0; i < nb_elements_f; i++)
+        {
+            const auto& element_indices = m_element_f_variable_indices[i];
+
+            std::vector<index> element_hi(length(element_indices));
+
+            for (index row_i = 0; row_i < length(element_indices); row_i++) {
+                const auto row = element_indices[row_i];
+                const auto col = element_indices.back();
+
+                element_hi[row_i] = m_structure_hm.get_index(row.global, col.global) + 1;
+            }
+
+            m_element_f_variable_indices_hi[i] = std::move(element_hi);
+        }
+
         Log::task_end("Problem initialized in {:.3f} sec", timer.ellapsed());
     }
 
@@ -457,16 +481,21 @@ class Problem {
 
             data.df(row.global) += g(row.local);
 
+            auto lo = m_structure_hm.get_first_index(row.global);
+            const auto hi = m_element_f_variable_indices_hi[i][row_i];
+
             for (index col_i = row_i; col_i < length(variable_indices) && TOrder > 1; col_i++) {
                 const auto col = variable_indices[col_i];
 
-                index index = m_structure_hm.get_index(row.global, col.global);
+                const index index = m_structure_hm.get_index_bounded(col.global, lo, hi);
 
                 if (row.local < col.local) {
                     data.hm_value(index) += h(row.local, col.local);
                 } else {
                     data.hm_value(index) += h(col.local, row.local);
                 }
+
+                lo = index;
             }
         }
 
@@ -597,9 +626,8 @@ class Problem {
             #pragma omp parallel if (m_nb_threads != 1) num_threads(m_nb_threads) firstprivate(l_data)
             {
                 #pragma omp for schedule(dynamic, m_grainsize) nowait
-                for (index i = 0; i < length(m_active_elements_f); i++) {
-                    // Log::info("Element {}", i);
-                    compute_element_f<TOrder>(l_data, m_active_elements_f[i]);
+                for (index i = 0; i < nb_elements_f(); i++) {
+                    compute_element_f<TOrder>(l_data, i);
                 }
 
                 if (sigma() != 1.0) {
@@ -615,8 +643,8 @@ class Problem {
                 }
 
                 #pragma omp for schedule(dynamic, m_grainsize) nowait
-                for (index i = 0; i < length(m_active_elements_g); i++) {
-                    compute_element_g<TOrder>(l_data, m_active_elements_g[i]);
+                for (index i = 0; i < nb_elements_g(); i++) {
+                    compute_element_g<TOrder>(l_data, i);
                 }
 
                 #pragma omp critical
diff --git a/include/eqlib/SparseStructure.h b/include/eqlib/SparseStructure.h
index acfd84f8..b42786ac 100644
--- a/include/eqlib/SparseStructure.h
+++ b/include/eqlib/SparseStructure.h
@@ -105,6 +105,31 @@ class SparseStructure {
         return m_ja.at(index);
     }
 
+    EQLIB_INLINE index get_first_index(const index i) const
+    {
+        return m_ia[i];
+    }
+
+    index get_index_bounded(const index j, const index lo, const index hi) const
+    {
+        const auto ja_begin = m_ja.begin();
+
+        const auto lower = std::next(ja_begin, lo);
+        const auto upper = std::next(ja_begin, hi);
+
+        const auto it = std::lower_bound(lower, upper, j);
+
+        if (*it != j || it == upper) {
+            return -1;
+        }
+
+        const index value_index = std::distance(ja_begin, it);
+
+        assert(value_index < nb_nonzeros());
+
+        return value_index;
+    }
+
     index get_index(const index row, const index col) const
     {
         assert(0 <= row && row < rows());
@@ -122,8 +147,10 @@ class SparseStructure {
 
             return it->second;
         } else {
-            const auto lower = m_ja.begin() + m_ia[i];
-            const auto upper = m_ja.begin() + m_ia[i + 1];
+            const auto ja_begin = m_ja.begin();
+
+            const auto lower = std::next(ja_begin, m_ia[i]);
+            const auto upper = std::next(ja_begin, m_ia[i + 1]);
 
             const auto it = std::lower_bound(lower, upper, j);
 
@@ -131,7 +158,7 @@ class SparseStructure {
                 return -1;
             }
 
-            const index value_index = std::distance(m_ja.begin(), it);
+            const index value_index = std::distance(ja_begin, it);
 
             assert(value_index < nb_nonzeros());