diff --git a/src/algebra/AlgebraicFunctions.cpp b/src/algebra/AlgebraicFunctions.cpp
index a2492884d10d276bc9db0196984b3196a7f64c97..36a75c69e860ccda4284b20a5649f2634273d099 100644
--- a/src/algebra/AlgebraicFunctions.cpp
+++ b/src/algebra/AlgebraicFunctions.cpp
@@ -209,12 +209,12 @@ namespace gmshfem::algebra
 
     if(A.format() == MatrixFormat::CRS) {
       const MatrixCRS< T_Scalar > &tmp = static_cast< const MatrixCRS< T_Scalar > & >(A);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
       {
         T_Scalar product = 0.;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < size; ++i) {
@@ -226,7 +226,7 @@ namespace gmshfem::algebra
           product += x[i] * line;
         }
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp critical
 #endif
         value += product;
@@ -234,12 +234,12 @@ namespace gmshfem::algebra
     }
     else if(A.format() == MatrixFormat::CRSFast) {
       const MatrixCRSFast< T_Scalar > &tmp = static_cast< const MatrixCRSFast< T_Scalar > & >(A);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
       {
         T_Scalar product = 0.;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < size; ++i) {
@@ -251,7 +251,7 @@ namespace gmshfem::algebra
           product += x[i] * line;
         }
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp critical
 #endif
         value += product;
@@ -279,11 +279,11 @@ namespace gmshfem::algebra
 
     if(A.format() == MatrixFormat::CRS) {
       const MatrixCRS< T_Scalar > &tmp = static_cast< const MatrixCRS< T_Scalar > & >(A);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
       {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < size; ++i) {
@@ -298,11 +298,11 @@ namespace gmshfem::algebra
     }
     else if(A.format() == MatrixFormat::CRSFast) {
       const MatrixCRSFast< T_Scalar > &tmp = static_cast< const MatrixCRSFast< T_Scalar > & >(A);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
       {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < size; ++i) {
@@ -335,12 +335,12 @@ namespace gmshfem::algebra
 
     if(A.format() == MatrixFormat::CRS) {
       const MatrixCRS< T_Scalar > &tmp = static_cast< const MatrixCRS< T_Scalar > & >(A);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
       {
         scalar::Precision< T_Scalar > product = 0.;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < size; ++i) {
@@ -352,7 +352,7 @@ namespace gmshfem::algebra
           product += std::abs((b[i] + line) * (b[i] + line));
         }
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp critical
 #endif
         value += product;
@@ -360,12 +360,12 @@ namespace gmshfem::algebra
     }
     else if(A.format() == MatrixFormat::CRSFast) {
       const MatrixCRSFast< T_Scalar > &tmp = static_cast< const MatrixCRSFast< T_Scalar > & >(A);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
       {
         scalar::Precision< T_Scalar > product = 0.;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < size; ++i) {
@@ -377,7 +377,7 @@ namespace gmshfem::algebra
           product += std::abs((b[i] + line) * (b[i] + line));
         }
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp critical
 #endif
         value += product;
diff --git a/src/algebra/MatrixCCS.cpp b/src/algebra/MatrixCCS.cpp
index 6b799de320da70d3f035e4b40197190136699e86..b8010f9491386f98c50d46f1bbea17208bc94348 100644
--- a/src/algebra/MatrixCCS.cpp
+++ b/src/algebra/MatrixCCS.cpp
@@ -87,11 +87,11 @@ namespace gmshfem::algebra
       this->copySize(other);
 
       _aj.resize(this->_size[1] + 1, 0);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
       for(auto i = 0ULL; i < mat.ai()[this->_size[0]]; ++i) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp atomic update
 #endif
         _aj[mat.aj()[i]]++;
diff --git a/src/algebra/MatrixCRS.cpp b/src/algebra/MatrixCRS.cpp
index 209fe20abfc1b3e98066f95626fecdf24d8d1ed4..366d1aff054a73b26a7cd2fb00b878b732f03da7 100644
--- a/src/algebra/MatrixCRS.cpp
+++ b/src/algebra/MatrixCRS.cpp
@@ -98,11 +98,11 @@ namespace gmshfem::algebra
       this->copySize(other);
 
       _ai.resize(this->_size[0] + 1, 0);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
       for(auto i = 0ULL; i < mat.aj()[this->_size[1]]; ++i) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp atomic update
 #endif
         _ai[mat.ai()[i]]++;
diff --git a/src/algebra/Vector.cpp b/src/algebra/Vector.cpp
index 889e676954b4ab6c73c6b576f130d644c1a59897..14e14e0079ac77d8506b74fb4379a9160279fb3f 100644
--- a/src/algebra/Vector.cpp
+++ b/src/algebra/Vector.cpp
@@ -278,7 +278,7 @@ namespace gmshfem::algebra
   {
     unsigned long long firstIndex = _vector.size();
     _vector.resize(_vector.size() + other._vector.size());
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
     for(auto i = 0ULL; i < other._vector.size(); ++i) {
diff --git a/src/analytics/Helmholtz2D.cpp b/src/analytics/Helmholtz2D.cpp
index 5bb445ad82ab8caa566656468cb7280466b4b5ad..d58a326993c33db0d99c47e8728427a21140c3b8 100644
--- a/src/analytics/Helmholtz2D.cpp
+++ b/src/analytics/Helmholtz2D.cpp
@@ -63,7 +63,7 @@ namespace gmshfem::analytics::helmholtz2D
       precomputed[n] = -(!n ? scalar::Precision< T_Scalar >(1.) : scalar::Precision< T_Scalar >(2.)) * JnkR0 * invHn1kR0 * powI[n % 4];
     }
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
     for(auto i = 0ULL; i < values.size(); ++i) {
@@ -125,7 +125,7 @@ namespace gmshfem::analytics::helmholtz2D
       precomputed[n] = -(!n ? scalar::Precision< T_Scalar >(1.) : scalar::Precision< T_Scalar >(2.)) * std::real(dHn1kR0) * invdHn1kR0 * powI[n % 4];
     }
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
     for(auto i = 0ULL; i < values.size(); ++i) {
@@ -187,7 +187,7 @@ namespace gmshfem::analytics::helmholtz2D
       precomputed[n] = -(!n ? scalar::Precision< T_Scalar >(1.) : scalar::Precision< T_Scalar >(2.)) * std::real(dHn1kR0) * invdHn1kR0 * powI[n % 4] * _k;
     }
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
     for(auto i = 0ULL; i < values.size(); ++i) {
@@ -250,7 +250,7 @@ namespace gmshfem::analytics::helmholtz2D
     else {
       kx = T_Scalar(-_M * _k / (1. - _M * _M), -std::sqrt(std::abs(ArgSqrt)) / (1. - _M * _M));
     }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
     for(auto i = 0ULL; i < values.size(); ++i) {
@@ -296,7 +296,7 @@ namespace gmshfem::analytics::helmholtz2D
     double valr, vali;
     const scalar::Precision< T_Scalar > ky = (_n * pi / _H);
     const scalar::Precision< T_Scalar > den = std::pow(_a * _k * _k, 2. / 3.);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
     for(auto i = 0ULL; i < values.size(); ++i) {
@@ -375,7 +375,7 @@ namespace gmshfem::analytics::helmholtz2D
       }
     }
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
     for(auto i = 0ULL; i < values.size(); ++i) {
@@ -431,7 +431,7 @@ namespace gmshfem::analytics::helmholtz2D
     const scalar::Precision< T_Scalar > Alphay = 1 + My * My / (beta * (1 + beta));
     const scalar::Precision< T_Scalar > K = Mx * My / (beta * (1 + beta));
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
     for(auto i = 0ULL; i < values.size(); ++i) {
diff --git a/src/analytics/Helmholtz3D.cpp b/src/analytics/Helmholtz3D.cpp
index 9abc430d4c834f461a64373d94adfe215f345457..3afbff28838a9a3984acd934a6c09b771c7d67b3 100644
--- a/src/analytics/Helmholtz3D.cpp
+++ b/src/analytics/Helmholtz3D.cpp
@@ -65,7 +65,7 @@ namespace gmshfem::analytics::helmholtz3D
     const scalar::Precision< T_Scalar > sinThetaInc = std::sin(_thetaInc);
     const scalar::Precision< T_Scalar > cosThetaInc = std::cos(_thetaInc);
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
     for(auto i = 0ULL; i < values.size(); ++i) {
@@ -141,7 +141,7 @@ namespace gmshfem::analytics::helmholtz3D
     const scalar::Precision< T_Scalar > sinThetaInc = std::sin(_thetaInc);
     const scalar::Precision< T_Scalar > cosThetaInc = std::cos(_thetaInc);
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
     for(auto i = 0ULL; i < values.size(); ++i) {
@@ -217,7 +217,7 @@ namespace gmshfem::analytics::helmholtz3D
     const scalar::Precision< T_Scalar > sinThetaInc = std::sin(_thetaInc);
     const scalar::Precision< T_Scalar > cosThetaInc = std::cos(_thetaInc);
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
     for(auto i = 0ULL; i < values.size(); ++i) {
@@ -294,7 +294,7 @@ namespace gmshfem::analytics::helmholtz3D
       precomputed[n] = powI[n % 4] * T_Scalar((2 * n + 1), 0.) * _kout * (djn1koutR0 * hn1koutR0 - dhn1koutR0 * jn1koutR0) / (rhokint * djn1kintR0 * hn1koutR0 - _kout * dhn1koutR0 * jn1kintR0);
     }
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
     for(auto i = 0ULL; i < values.size(); ++i) {
diff --git a/src/analytics/Maxwell2D.cpp b/src/analytics/Maxwell2D.cpp
index 8f268dd844b886579fd0b60f2f4cdb78fb227989..a1e316bf877b44dd18125ecc903a1662842d77a9 100644
--- a/src/analytics/Maxwell2D.cpp
+++ b/src/analytics/Maxwell2D.cpp
@@ -48,14 +48,14 @@ namespace gmshfem::analytics::maxwell2D
   template< class T_Scalar >
   void ScatteringByAPECCylinderTE< T_Scalar >::operator()(function::OutputVector< T_Scalar, Degree::Degree1 > &values, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points) const
   {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
     _values.resize(values.size());
 
     _helmholtz(_values, points);
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
     for(auto i = 0ULL; i < values.size(); ++i) {
@@ -108,7 +108,7 @@ namespace gmshfem::analytics::maxwell2D
       precomputed[n] = (!n ? scalar::Precision< T_Scalar >(1.) : scalar::Precision< T_Scalar >(2.)) * std::real(dHn1kR0) * invdHn1kR0 * (!n ? powI[3] : powI[(n - 1) % 4]) / _k;
     }
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
     for(auto i = 0ULL; i < values.size(); ++i) {
@@ -166,14 +166,14 @@ namespace gmshfem::analytics::maxwell2D
   template< class T_Scalar >
   void ScatteringByAPMCCylinderTE< T_Scalar >::operator()(function::OutputVector< T_Scalar, Degree::Degree1 > &values, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points) const
   {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
     _values.resize(values.size());
 
     _helmholtz(_values, points);
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
     for(auto i = 0ULL; i < values.size(); ++i) {
@@ -216,14 +216,14 @@ namespace gmshfem::analytics::maxwell2D
   template< class T_Scalar >
   void dr_ScatteringByAPMCCylinderTE< T_Scalar >::operator()(function::OutputVector< T_Scalar, Degree::Degree1 > &values, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points) const
   {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
     _values.resize(values.size());
 
     _helmholtz(_values, points);
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
     for(auto i = 0ULL; i < values.size(); ++i) {
diff --git a/src/analytics/Maxwell3D.cpp b/src/analytics/Maxwell3D.cpp
index 2577e2ba1fe2ab68e4294b2c2961fa7bed7c7761..3dd03888694a41ab93b1d61efeaa46e2f166dee3 100644
--- a/src/analytics/Maxwell3D.cpp
+++ b/src/analytics/Maxwell3D.cpp
@@ -69,7 +69,7 @@ namespace gmshfem::analytics::maxwell3D
       precomputedBlm[n - 1] = powI[(n + 1) % 4] * T_Scalar((2. * n + 1) / (n * (n + 1)), 0.) * (_k2_out * _k_int * koutR0jn1koutR0 * dkintR0jn1kintR0 - _k2_int * _k_out * dkoutR0jn1koutR0 * kintR0jn1kintR0) / (_k2_out * _k_int * koutR0hn1koutR0 * dkintR0jn1kintR0 - _k2_int * _k_out * dkoutR0hn1koutR0 * kintR0jn1kintR0);
     }
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
     for(auto i = 0ULL; i < values.size(); ++i) {
@@ -174,7 +174,7 @@ namespace gmshfem::analytics::maxwell3D
       precomputedBlm[n - 1] = powI[(n + 1) % 4] * T_Scalar((2. * n + 1) / (n * (n + 1)), 0.) * (_k2_out * _k_int * koutR0jn1koutR0 * dkintR0jn1kintR0 - _k2_int * _k_out * dkoutR0jn1koutR0 * kintR0jn1kintR0) / (_k2_out * _k_int * koutR0hn1koutR0 * dkintR0jn1kintR0 - _k2_int * _k_out * dkoutR0hn1koutR0 * kintR0jn1kintR0);
     }
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
     for(auto i = 0ULL; i < values.size(); ++i) {
diff --git a/src/analytics/Navier2D.cpp b/src/analytics/Navier2D.cpp
index 4ec4792cdea8fda63e9cb9eb701361fcd9e745fe..bacbb0c51ba92221cdc39592ebf91391f5b406ae 100644
--- a/src/analytics/Navier2D.cpp
+++ b/src/analytics/Navier2D.cpp
@@ -126,7 +126,7 @@ namespace gmshfem::analytics::navier2D
       B[n] = scalar::Precision< T_Scalar >(1.) / (m11 * m22 - m12 * m21) * (-m21 * f_1 + m11 * f_2);
     }
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
     for(auto i = 0ULL; i < values.size(); ++i) {
@@ -202,7 +202,7 @@ namespace gmshfem::analytics::navier2D
       B[n] = scalar::Precision< T_Scalar >(1.) / (m11 * m22 - m12 * m21) * (-m21 * f_1 + m11 * f_2);
     }
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
     for(auto i = 0ULL; i < values.size(); ++i) {
@@ -327,7 +327,7 @@ namespace gmshfem::analytics::navier2D
       AB4[n] = inverseM41 * f1 + inverseM42 * f2 + inverseM43 * f3 + inverseM44 * f4;
     }
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
     for(auto i = 0ULL; i < values.size(); ++i) {
@@ -450,7 +450,7 @@ namespace gmshfem::analytics::navier2D
       AB4[n] = inverseM41 * f1 + inverseM42 * f2 + inverseM43 * f3 + inverseM44 * f4;
     }
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
     for(auto i = 0ULL; i < values.size(); ++i) {
@@ -610,7 +610,7 @@ namespace gmshfem::analytics::navier2D
       AB4[n] = inverseM41 * f1 + inverseM42 * f2 + inverseM43 * f3 + inverseM44 * f4;
     }
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
     for(auto i = 0ULL; i < values.size(); ++i) {
@@ -770,7 +770,7 @@ namespace gmshfem::analytics::navier2D
       AB4[n] = inverseM41 * f1 + inverseM42 * f2 + inverseM43 * f3 + inverseM44 * f4;
     }
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
     for(auto i = 0ULL; i < values.size(); ++i) {
diff --git a/src/common/io/PPMio.cpp b/src/common/io/PPMio.cpp
index 0824986151d482060d447c5aa58435827e66ec02..eb287793d571109448bce34485dc83404c1c6095 100644
--- a/src/common/io/PPMio.cpp
+++ b/src/common/io/PPMio.cpp
@@ -115,7 +115,7 @@ namespace gmshfem::common
       if(_tmp.size() != 3 * map.size()) {
         _tmp.resize(3 * map.size());
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
       for(auto i = 0ULL; i < map.size(); ++i) {
@@ -126,7 +126,7 @@ namespace gmshfem::common
       if(_tmp.size() != map.size()) {
         _tmp.resize(map.size());
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
       for(auto i = 0ULL; i < map.size(); ++i) {
@@ -137,13 +137,13 @@ namespace gmshfem::common
       if(_tmp.size() != (map.size() + 7) / 8) {
         _tmp.resize((map.size() + 7) / 8);
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
       for(auto i = 0ULL; i < _tmp.size(); ++i) {
         _tmp[i] = 0;
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
       for(auto i = 0ULL; i < map.size(); ++i) {
diff --git a/src/common/numa.h b/src/common/numa.h
index 1f5dc4c4d6dda5873c7bb08d72de7b26b3b03851..6a6752b46f489ebb6f5723c539c45a84fb920e35 100644
--- a/src/common/numa.h
+++ b/src/common/numa.h
@@ -38,7 +38,7 @@ namespace gmshfem::numa
       T *p = std::allocator< T >::allocate(n);
 
       // first touch
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
       for(unsigned long long i = 0; i < n; ++i) {
@@ -58,7 +58,7 @@ namespace gmshfem::numa
   void copy(std::vector< T_Scalar1, allocator< T_Scalar1 > > &v1, std::vector< T_Scalar2 > &v2)
   {
     v1.resize(v2.size());
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
     for(auto i = 0ULL; i < v2.size(); ++i) {
diff --git a/src/common/scalar.h b/src/common/scalar.h
index e644ac774f8160a9261162fdc72b251e8e0e9feb..48ee0985d3c446778f6dd1ccf3734c4a75e88b3e 100644
--- a/src/common/scalar.h
+++ b/src/common/scalar.h
@@ -22,6 +22,9 @@ namespace Eigen
 #include <limits>
 #include <string>
 #include <vector>
+#ifdef HAVE_PETSC
+#include <petsc.h>
+#endif // HAVE_PETSC
 #ifdef DDMGPU
 #include <thrust/complex.h>
 #endif
diff --git a/src/dofs/DofsFactory.cpp b/src/dofs/DofsFactory.cpp
index 5b2ae014f167bc65db40392d8ce244d65e998afc..ad4a5f39d692dfb3922c7e4c4aa432800bdb382c 100644
--- a/src/dofs/DofsFactory.cpp
+++ b/src/dofs/DofsFactory.cpp
@@ -81,7 +81,7 @@ namespace gmshfem::dofs
               std::vector< scalar::Precision< T_Scalar > > coordNoNuma;
               field->getFunctionSpace()->getKeys(true, typeKeys, entityKeys, coordNoNuma, elementTypes[typeIndex], *itEntity);
               numa::copy(coord, coordNoNuma);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
               {
@@ -188,7 +188,7 @@ namespace gmshfem::dofs
               std::vector< scalar::Precision< T_Scalar > > coordNoNuma;
               field->getFunctionSpace()->getKeys(true, typeKeys, entityKeys, coordNoNuma, elementTypes[typeIndex], *itEntity);
               numa::copy(coord, coordNoNuma);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
   #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
               {
@@ -295,7 +295,7 @@ namespace gmshfem::dofs
               std::vector< scalar::Precision< T_Scalar > > coordNoNuma;
               field->getFunctionSpace()->getKeys(true, typeKeys, entityKeys, coordNoNuma, elementTypes[typeIndex], *itEntity);
               numa::copy(coord, coordNoNuma);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
   #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
               {
diff --git a/src/dofs/DofsManager.cpp b/src/dofs/DofsManager.cpp
index d76af7519c98b5fdee57e161686b7511455b9e80..9d3a9e9651a69f3090ec00a8944713aca040426e 100644
--- a/src/dofs/DofsManager.cpp
+++ b/src/dofs/DofsManager.cpp
@@ -177,7 +177,7 @@ namespace gmshfem::dofs
             static_cast< UnknownDof * >(sortedEntitiesUnknownDofs[i].ptrB)->numDof(i + 1 - globalDofShift);
           }
         }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
         for(auto i = 0ULL; i < sortedEntitiesUnknownBubbleDofs.size(); ++i) {
@@ -236,7 +236,7 @@ namespace gmshfem::dofs
             static_cast< UnknownDof * >(sortedEntitiesUnknownDofs[i].ptrB)->numDof(i + 1 - globalDofShift);
           }
         }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
         for(auto i = 0ULL; i < sortedEntitiesUnknownBubbleDofs.size(); ++i) {
diff --git a/src/domain/Axisymmetry.cpp b/src/domain/Axisymmetry.cpp
index 8f9fb0a5f831fe95e29c238d5408c25994690bdb..4d1b37a16d253a16725668628b79d4f1f9f983f4 100644
--- a/src/domain/Axisymmetry.cpp
+++ b/src/domain/Axisymmetry.cpp
@@ -35,7 +35,7 @@ namespace gmshfem::domain
   {
     const unsigned long long nbrValues = points.size() / 3;
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
     for(auto i = 0ULL; i < nbrValues; ++i) {
@@ -43,7 +43,7 @@ namespace gmshfem::domain
 
       if(R == 0.) {
         R = scalar::Epsilon< T_PScalar >::value;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp critical
 #endif
         msg::debug << "A zero axisymmetric shell coordinate is shifted by eps0 (" << R << ")" << msg::endl;
diff --git a/src/domain/AxisymmetryShell.cpp b/src/domain/AxisymmetryShell.cpp
index ec82948358a419a3f93a654f7dc36d43d1119b18..2d3bf7422673524b0df3269f0223a16b7338dbd7 100644
--- a/src/domain/AxisymmetryShell.cpp
+++ b/src/domain/AxisymmetryShell.cpp
@@ -37,7 +37,7 @@ namespace gmshfem::domain
   {
     const unsigned long long nbrValues = points.size() / 3;
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
     for(auto i = 0ULL; i < nbrValues; ++i) {
@@ -54,7 +54,7 @@ namespace gmshfem::domain
 
       if(S == 0.) {
         S = scalar::Epsilon< T_PScalar >::value;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp critical
 #endif
         msg::debug << "A zero axisymmetric shell coordinate is shifted by eps0(" << S << ")" << msg::endl;
diff --git a/src/domain/CylindricalShell.cpp b/src/domain/CylindricalShell.cpp
index bb75c849970f81e4e60cff86a46cd8cef2f643f8..99cda7b7e10eab0afc4966c61db0dae6fba8b64b 100644
--- a/src/domain/CylindricalShell.cpp
+++ b/src/domain/CylindricalShell.cpp
@@ -39,7 +39,7 @@ namespace gmshfem::domain
   {
     const unsigned long long nbrValues = points.size() / 3;
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
     for(auto i = 0ULL; i < nbrValues; ++i) {
diff --git a/src/domain/Domain.cpp b/src/domain/Domain.cpp
index bfbd485e2f2eabaeec88758b0dfdd9122acbf8b3..06d9a491667170a3f9c01b47b11335b354e79876 100644
--- a/src/domain/Domain.cpp
+++ b/src/domain/Domain.cpp
@@ -301,27 +301,27 @@ namespace gmshfem::domain
     for(auto it = cbegin(); it != cend(); ++it) {
       gmsh::model::mesh::getElementTypes(elementTypes, it->first, it->second);
       for(unsigned int iType = 0; iType < elementTypes.size(); ++iType) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
         {
           const unsigned int numThreads = omp::getNumThreads();
           const unsigned int myThreadID = omp::getThreadNum();
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
           gmsh::model::mesh::preallocateElementsByType(elementTypes[iType], true, true, elementTags, nodeTags, it->second);
           gmsh::model::mesh::getElementsByType(elementTypes[iType], elementTags, nodeTags, it->second, myThreadID, numThreads);
           const unsigned int nodeByElement = nodeTags.size() / elementTags.size();
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
           for(unsigned int i = 0; i < elementTags.size(); ++i) {
             for(unsigned int j = 0; j < nodeByElement; ++j) {
               auto itFind = toolNodes.find(nodeTags[nodeByElement * i + j]);
               if(itFind != toolNodes.end()) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp critical
 #endif
                 {
diff --git a/src/domain/PolarShell.cpp b/src/domain/PolarShell.cpp
index a0305ab0410480d047b3ff42382d1f67ec26aebe..1b099c2824c881788e9bca7f1333614778c01584 100644
--- a/src/domain/PolarShell.cpp
+++ b/src/domain/PolarShell.cpp
@@ -39,7 +39,7 @@ namespace gmshfem::domain
   {
     const unsigned long long nbrValues = points.size() / 3;
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
     for(auto i = 0ULL; i < nbrValues; ++i) {
diff --git a/src/domain/SphericalShell.cpp b/src/domain/SphericalShell.cpp
index 5129a6bac6271f0364cfdcfd6b95e9f82bb78691..f3dcfac58a53e0c18acc7757fe375da9e85732c2 100644
--- a/src/domain/SphericalShell.cpp
+++ b/src/domain/SphericalShell.cpp
@@ -39,7 +39,7 @@ namespace gmshfem::domain
   {
     const unsigned long long nbrValues = points.size() / 3;
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
     for(auto i = 0ULL; i < nbrValues; ++i) {
diff --git a/src/field/FieldInterface.cpp b/src/field/FieldInterface.cpp
index 32d15e572d2184abe4a68ee558ef1c8b6f26cfde..87d48e28a03a5d3de4414a705ce7243c6159105e 100644
--- a/src/field/FieldInterface.cpp
+++ b/src/field/FieldInterface.cpp
@@ -1091,7 +1091,7 @@ namespace gmshfem::field
         std::vector< unsigned long long > numIndices(it->second.size());
         std::vector< int > typeIndices(it->second.size());
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
         for(auto i = 0ULL; i < it->second.size(); ++i) {
@@ -1127,7 +1127,7 @@ namespace gmshfem::field
         }
         dofs.resize(typeKeys.size() * multiplicity());
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
         {
@@ -1156,7 +1156,7 @@ namespace gmshfem::field
     if(!indices.have(_tag)) {
       std::vector< unsigned long long > numIndices(multiplicity() * typeKeys.size());
       std::vector< int > typeIndices(multiplicity() * typeKeys.size());
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
       {
diff --git a/src/field/FunctionSpaceConstant.cpp b/src/field/FunctionSpaceConstant.cpp
index 3d6f8d6df4364659afd3a0dc79ede31bda19862d..c4ed7577448982d4b757a76f9e68475c99f15e5b 100644
--- a/src/field/FunctionSpaceConstant.cpp
+++ b/src/field/FunctionSpaceConstant.cpp
@@ -120,11 +120,11 @@ namespace gmshfem::field
     typeKeys.resize(elements.size());
     entityKeys.resize(elements.size());
     std::vector< double > gmshCoord;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < elements.size(); ++i) {
@@ -136,13 +136,13 @@ namespace gmshfem::field
         const unsigned int numThreads = omp::getNumThreads();
         const unsigned int myThreadID = omp::getThreadNum();
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         gmsh::model::mesh::preallocateBarycenters(elementType, gmshCoord, entity.second);
 
         gmsh::model::mesh::getBarycenters(elementType, entity.second, false, true, gmshCoord, myThreadID, numThreads);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #pragma omp single
 #endif
diff --git a/src/function/ScalarFunction.cpp b/src/function/ScalarFunction.cpp
index dbb7cd2ed4aae7139497a07e55d51688a78dbf9c..18ea0ad03482af901007b4b0613bef6abbb8406f 100644
--- a/src/function/ScalarFunction.cpp
+++ b/src/function/ScalarFunction.cpp
@@ -207,7 +207,7 @@ namespace gmshfem::function
     }
     if(_tree) {
       if(values.size() < points.size() / 3) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #pragma omp single
 #endif
@@ -215,11 +215,11 @@ namespace gmshfem::function
       }
       FunctionAllocator< typename MathObject< T_Scalar, Degree::Degree0 >::Object > allocator(&values[0]);
       OutputVector< T_Scalar, Degree::Degree0 > valuesOV(points.size() / 3, allocator);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
       _tree->evaluate(valuesOV, points, gaussPoints, elementType, entity);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 
 #pragma omp single
@@ -266,7 +266,7 @@ namespace gmshfem::function
 
           gmsh::model::mesh::preallocateJacobians(elementTypes[typeIndex], nbrOfGaussPoints, false, true, needPoints, gmshJacobians, gmshDeterminants, gmshPoints, itEntity->second);
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
           {
@@ -289,7 +289,7 @@ namespace gmshfem::function
           FunctionAllocator< typename MathObject< T_Scalar, Degree::Degree0 >::Object > allocator(&values[i][0]);
           OutputVector< T_Scalar, Degree::Degree0 > valuesOV(points.size() / 3, allocator);
           timer.tick();
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
           _tree->evaluate(valuesOV, points, gaussPoints, elementTypes[typeIndex], *itEntity);
diff --git a/src/function/Tensor4Function.cpp b/src/function/Tensor4Function.cpp
index 8778d53c38aa13f5da1330b8217cc5fba7f69480..3cfe8837d386fffed6fc057a6f5031756be65a80 100644
--- a/src/function/Tensor4Function.cpp
+++ b/src/function/Tensor4Function.cpp
@@ -155,7 +155,7 @@ namespace gmshfem::function
     }
     if(_tree) {
       if(values.size() < points.size() / 3) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 
 #pragma omp single
@@ -164,11 +164,11 @@ namespace gmshfem::function
       }
       FunctionAllocator< typename MathObject< T_Scalar, Degree::Degree4 >::Object > allocator(&values[0]);
       OutputVector< T_Scalar, Degree::Degree4 > valuesOV(points.size() / 3, allocator);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
       _tree->evaluate(valuesOV, points, gaussPoints, elementType, entity);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 
 #pragma omp single
@@ -215,7 +215,7 @@ namespace gmshfem::function
 
           gmsh::model::mesh::preallocateJacobians(elementTypes[typeIndex], nbrOfGaussPoints, false, true, needPoints, gmshJacobians, gmshDeterminants, gmshPoints, itEntity->second);
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
           {
@@ -238,7 +238,7 @@ namespace gmshfem::function
           FunctionAllocator< typename MathObject< T_Scalar, Degree::Degree4 >::Object > allocator(&values[i][0]);
           OutputVector< T_Scalar, Degree::Degree4 > valuesOV(points.size() / 3, allocator);
           timer.tick();
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
           _tree->evaluate(valuesOV, points, gaussPoints, elementTypes[typeIndex], *itEntity);
diff --git a/src/function/TensorFunction.cpp b/src/function/TensorFunction.cpp
index 109dbb774ac0dcd2ab1732334b3d95471a0629ac..6286010c5ff8a8c0604a9271fe0c5d2b9d98ebbe 100644
--- a/src/function/TensorFunction.cpp
+++ b/src/function/TensorFunction.cpp
@@ -223,7 +223,7 @@ namespace gmshfem::function
     }
     if(_tree) {
       if(values.size() < points.size() / 3) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 
 #pragma omp single
@@ -232,11 +232,11 @@ namespace gmshfem::function
       }
       FunctionAllocator< typename MathObject< T_Scalar, Degree::Degree2 >::Object > allocator(&values[0]);
       OutputVector< T_Scalar, Degree::Degree2 > valuesOV(points.size() / 3, allocator);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
       _tree->evaluate(valuesOV, points, gaussPoints, elementType, entity);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 
 #pragma omp single
@@ -283,7 +283,7 @@ namespace gmshfem::function
 
           gmsh::model::mesh::preallocateJacobians(elementTypes[typeIndex], nbrOfGaussPoints, false, true, needPoints, gmshJacobians, gmshDeterminants, gmshPoints, itEntity->second);
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
           {
@@ -306,7 +306,7 @@ namespace gmshfem::function
           FunctionAllocator< typename MathObject< T_Scalar, Degree::Degree2 >::Object > allocator(&values[i][0]);
           OutputVector< T_Scalar, Degree::Degree2 > valuesOV(points.size() / 3, allocator);
           timer.tick();
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
           _tree->evaluate(valuesOV, points, gaussPoints, elementTypes[typeIndex], *itEntity);
diff --git a/src/function/VectorFunction.cpp b/src/function/VectorFunction.cpp
index 32de97f85f810da0ea550e89fb34c45e5ccb835c..ebbe484c901feabd7fb853070ce48507a5a39072 100644
--- a/src/function/VectorFunction.cpp
+++ b/src/function/VectorFunction.cpp
@@ -257,7 +257,7 @@ namespace gmshfem::function
     }
     if(_tree) {
       if(values.size() < points.size() / 3) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 
 #pragma omp single
@@ -266,11 +266,11 @@ namespace gmshfem::function
       }
       FunctionAllocator< typename MathObject< T_Scalar, Degree::Degree1 >::Object > allocator(&values[0]);
       OutputVector< T_Scalar, Degree::Degree1 > valuesOV(points.size() / 3, allocator);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
       _tree->evaluate(valuesOV, points, gaussPoints, elementType, entity);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 
 #pragma omp single
@@ -317,7 +317,7 @@ namespace gmshfem::function
 
           gmsh::model::mesh::preallocateJacobians(elementTypes[typeIndex], nbrOfGaussPoints, false, true, needPoints, gmshJacobians, gmshDeterminants, gmshPoints, itEntity->second);
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
           {
@@ -340,7 +340,7 @@ namespace gmshfem::function
           FunctionAllocator< typename MathObject< T_Scalar, Degree::Degree1 >::Object > allocator(&values[i][0]);
           OutputVector< T_Scalar, Degree::Degree1 > valuesOV(points.size() / 3, allocator);
           timer.tick();
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
           _tree->evaluate(valuesOV, points, gaussPoints, elementTypes[typeIndex], *itEntity);
diff --git a/src/function/executionTree/AnalyticalNode.h b/src/function/executionTree/AnalyticalNode.h
index 7fdaf61f78928c209070d58caa6e4d22bed212bd..873b684edcc9c9681083034d35bc56eb10bdb7b6 100644
--- a/src/function/executionTree/AnalyticalNode.h
+++ b/src/function/executionTree/AnalyticalNode.h
@@ -79,20 +79,20 @@ namespace gmshfem::function
     virtual void evaluate(OutputVector< typename T_Op::D_Scalar, T_Op::D_Degree > &values, const std::vector< scalar::Precision< typename T_Op::D_Scalar >, numa::allocator< scalar::Precision< typename T_Op::D_Scalar > > > &points, const std::vector< scalar::Precision< typename T_Op::D_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const override
     {
       if(values.size() < points.size() / 3) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
         FunctionAllocator< typename MathObject< typename T_Op::D_Scalar, T_Op::D_Degree >::Object >::initalizeNUMAMemory(points.size() / 3);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp master
 #endif
         values.resize(points.size() / 3);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
       }
       _op(values, points);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
     }
diff --git a/src/function/executionTree/BinaryNode.h b/src/function/executionTree/BinaryNode.h
index b58f17db541978d06fd2d3f534565288e7542e47..5579f8e6f9fdd8ec4793d60c5e9f119a863701bd 100644
--- a/src/function/executionTree/BinaryNode.h
+++ b/src/function/executionTree/BinaryNode.h
@@ -101,7 +101,7 @@ namespace gmshfem::function
       {
         bool isAllConstant = true;
         for(int i = 0; i < 2; ++i) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp critical
 #endif
           enclose._isConstant[i] = enclose._leaves.get(i)->isConstant(entity);
@@ -170,7 +170,7 @@ namespace gmshfem::function
             }
             return;
           });
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
         common::turn_into_tuple< 0, 2 >(
@@ -299,21 +299,21 @@ namespace gmshfem::function
     virtual void evaluate(OutputVector< typename T_Op::D_Scalar, T_Op::D_Degree > &values, const std::vector< scalar::Precision< typename T_Op::D_Scalar >, numa::allocator< scalar::Precision< typename T_Op::D_Scalar > > > &points, const std::vector< scalar::Precision< typename T_Op::D_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const override
     {
       if(values.size() < points.size() / 3) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
         FunctionAllocator< typename MathObject< typename T_Op::D_Scalar, T_Op::D_Degree >::Object >::initalizeNUMAMemory(points.size() / 3);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp master
 #endif
         values.resize(points.size() / 3);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
       }
       EfficientEvaluate evaluator;
       evaluator(*this, values, points, gaussPoints, elementType, entity);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 
 #pragma omp single
@@ -322,7 +322,7 @@ namespace gmshfem::function
         _tmpA.clear();
         _tmpA.shrink_to_fit();
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
diff --git a/src/function/executionTree/ChangeOfCoordinatesNode.h b/src/function/executionTree/ChangeOfCoordinatesNode.h
index 81b54b25b53cf0a065c3395a9218372ba4c372a7..5f5a99e729dd5c7cada663580e001a61f9a975e6 100644
--- a/src/function/executionTree/ChangeOfCoordinatesNode.h
+++ b/src/function/executionTree/ChangeOfCoordinatesNode.h
@@ -107,22 +107,22 @@ namespace gmshfem::function
 
     virtual void evaluate(OutputVector< T_Scalar, T_Degree > &values, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const override
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp master
 #endif
       _tmpPoints.resize(points.size());
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
 
       for(auto i = 0; i < 3; ++i) {
         _coordLeaf[i]->evaluate(_tmpCoord, _coordLeaf[i]->isConstant(entity) ? _zeros : points, gaussPoints, elementType, entity);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
         if(_coordLeaf[i]->isConstant(entity)) {
           InputVector< scalar::Precision< T_Scalar >, Degree::Degree0, std::true_type > inputVector(_tmpCoord);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
           for(auto j = 0ULL; j < points.size() / 3; ++j) {
@@ -131,30 +131,30 @@ namespace gmshfem::function
         }
         else {
           InputVector< scalar::Precision< T_Scalar >, Degree::Degree0, std::false_type > inputVector(_tmpCoord);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
           for(auto j = 0ULL; j < points.size() / 3; ++j) {
             _tmpPoints[3 * j + i] = inputVector[j];
           }
         }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         _tmpCoord.clear();
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
         _tmpCoord.clear();
         _tmpCoord.shrink_to_fit();
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp master
 #endif
       values.resize(points.size() / 3);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
 
diff --git a/src/function/executionTree/FieldNode.h b/src/function/executionTree/FieldNode.h
index f3e76e8fcc81271bf8a2fa70e7785b176d1d2638..453b22a7b21b123a5015ec75bb8921be30847b45 100644
--- a/src/function/executionTree/FieldNode.h
+++ b/src/function/executionTree/FieldNode.h
@@ -81,21 +81,21 @@ namespace gmshfem::function
     virtual void evaluate(OutputVector< typename T_Op::D_Scalar, T_Op::D_Degree > &values, const std::vector< scalar::Precision< typename T_Op::D_Scalar >, numa::allocator< scalar::Precision< typename T_Op::D_Scalar > > > &points, const std::vector< scalar::Precision< typename T_Op::D_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const override
     {
       if(values.size() < points.size() / 3) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
         FunctionAllocator< typename MathObject< typename T_Op::D_Scalar, T_Op::D_Degree >::Object >::initalizeNUMAMemory(points.size() / 3);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp master
 #endif
         values.resize(points.size() / 3);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
       }
 
       _op(_field, values, points, gaussPoints, elementType, entity);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
     }
@@ -206,21 +206,21 @@ namespace gmshfem::function
     virtual void evaluate(OutputVector< typename T_Op::D_ScalarOut, T_Op::D_Degree > &values, const std::vector< scalar::Precision< typename T_Op::D_ScalarIn >, numa::allocator< scalar::Precision< typename T_Op::D_ScalarIn > > > &points, const std::vector< scalar::Precision< typename T_Op::D_ScalarIn > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const override
     {
       if(values.size() < points.size() / 3) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
         FunctionAllocator< typename MathObject< typename T_Op::D_ScalarOut, T_Op::D_Degree >::Object >::initalizeNUMAMemory(points.size() / 3);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp master
 #endif
         values.resize(points.size() / 3);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
       }
 
       _op(_field, values, points, gaussPoints, elementType, entity);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
     }
diff --git a/src/function/executionTree/FunctionAllocator.h b/src/function/executionTree/FunctionAllocator.h
index 118f285364b2d039ac3b782ec2e754f13ecaa84d..f3db33ee14cb4c220d0b27683a26e705b916c634 100644
--- a/src/function/executionTree/FunctionAllocator.h
+++ b/src/function/executionTree/FunctionAllocator.h
@@ -140,26 +140,26 @@ namespace gmshfem::function
       }
       else {
         static T *s_p = nullptr;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp master
 #endif
         s_p = static_cast< T * >(MemoryPoolAllocator::instance()->allocate(n * sizeof(value_type)));
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(size_type i = 0; i < n; ++i) {
           *reinterpret_cast< char * >(s_p + i) = 0.;
         }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 
 #pragma omp master
 #endif
         MemoryPoolAllocator::instance()->deallocate(s_p, n * sizeof(value_type));
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
       }
diff --git a/src/function/executionTree/MultaryNode.h b/src/function/executionTree/MultaryNode.h
index d1c0eda36446e777edfb287c8c3ff69ed518c96f..370a349b20c6467169c2b44922d83a82264fb065 100644
--- a/src/function/executionTree/MultaryNode.h
+++ b/src/function/executionTree/MultaryNode.h
@@ -60,19 +60,19 @@ namespace gmshfem::function
     virtual void evaluate(OutputVector< T_Scalar, T_Degree > &values, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const override
     {
       if(values.size() < points.size() / 3) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
         FunctionAllocator< typename MathObject< T_Scalar, T_Degree >::Object >::initalizeNUMAMemory(points.size() / 3);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp master
 #endif
         values.resize(points.size() / 3);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
@@ -85,7 +85,7 @@ namespace gmshfem::function
       if(_it != _leaves.end()) {
         _it->second->evaluate(values, points, gaussPoints, elementType, entity);
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
     }
diff --git a/src/function/executionTree/NovenaryNode.h b/src/function/executionTree/NovenaryNode.h
index e0aab867e48f7bcdbadce1985e20e62fc869a451..295c3664995a486bd3adceea0dfaa6385c850dab 100644
--- a/src/function/executionTree/NovenaryNode.h
+++ b/src/function/executionTree/NovenaryNode.h
@@ -158,7 +158,7 @@ namespace gmshfem::function
       {
         bool isAllConstant = true;
         for(int i = 0; i < 9; ++i) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp critical
 #endif
           enclose._isConstant[i] = enclose._leaves.get(i)->isConstant(entity);
@@ -265,7 +265,7 @@ namespace gmshfem::function
             return;
           });
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
         common::turn_into_tuple< 0, 9 >(
@@ -404,21 +404,21 @@ namespace gmshfem::function
     virtual void evaluate(OutputVector< typename T_Op::D_Scalar, T_Op::D_Degree > &values, const std::vector< scalar::Precision< typename T_Op::D_Scalar >, numa::allocator< scalar::Precision< typename T_Op::D_Scalar > > > &points, const std::vector< scalar::Precision< typename T_Op::D_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const override
     {
       if(values.size() < points.size() / 3) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
         FunctionAllocator< typename MathObject< typename T_Op::D_Scalar, T_Op::D_Degree >::Object >::initalizeNUMAMemory(points.size() / 3);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp master
 #endif
         values.resize(points.size() / 3);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
       }
       EfficientEvaluate evaluator;
       evaluator(*this, values, points, gaussPoints, elementType, entity);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 
 #pragma omp single
@@ -427,56 +427,56 @@ namespace gmshfem::function
         _tmpA.clear();
         _tmpA.shrink_to_fit();
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
         _tmpB.clear();
         _tmpB.shrink_to_fit();
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
         _tmpC.clear();
         _tmpC.shrink_to_fit();
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
         _tmpD.clear();
         _tmpD.shrink_to_fit();
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
         _tmpE.clear();
         _tmpE.shrink_to_fit();
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
         _tmpF.clear();
         _tmpF.shrink_to_fit();
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
         _tmpG.clear();
         _tmpG.shrink_to_fit();
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
         _tmpH.clear();
         _tmpH.shrink_to_fit();
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
diff --git a/src/function/executionTree/NullaryNode.h b/src/function/executionTree/NullaryNode.h
index 15c395866710f68f82563567660627db2ba21d55..69bdfb14da9f75e387f7885e9cdbffc27397f776 100644
--- a/src/function/executionTree/NullaryNode.h
+++ b/src/function/executionTree/NullaryNode.h
@@ -76,20 +76,20 @@ namespace gmshfem::function
     virtual void evaluate(OutputVector< typename T_Op::D_Scalar, T_Op::D_Degree > &values, const std::vector< scalar::Precision< typename T_Op::D_Scalar >, numa::allocator< scalar::Precision< typename T_Op::D_Scalar > > > &points, const std::vector< scalar::Precision< typename T_Op::D_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const override
     {
       if(values.size() < points.size() / 3) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
         FunctionAllocator< typename MathObject< typename T_Op::D_Scalar, T_Op::D_Degree >::Object >::initalizeNUMAMemory(points.size() / 3);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp master
 #endif
         values.resize(points.size() / 3);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
       }
       _op(values, points, gaussPoints, elementType, entity);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
     }
diff --git a/src/function/executionTree/ScalarTypeNode.h b/src/function/executionTree/ScalarTypeNode.h
index b84d6902f799d1531e88392da759c63f8750e0ab..daa89e074a10ae55edc4f5ae9c586d2149472b97 100644
--- a/src/function/executionTree/ScalarTypeNode.h
+++ b/src/function/executionTree/ScalarTypeNode.h
@@ -77,21 +77,21 @@ namespace gmshfem::function
     {
       this->_leaf.get(0)->evaluate(_tmpA, points, gaussPoints, elementType, entity);
       if(values.size() < points.size() / 3) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
         FunctionAllocator< typename MathObject< typename T_Op::D_ScalarOut, T_Op::D_Degree >::Object >::initalizeNUMAMemory(points.size() / 3);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp master
 #endif
         values.resize(points.size() / 3);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
       }
 
       _op(values, _tmpA, points, gaussPoints, elementType, entity);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 
 #pragma omp single
diff --git a/src/function/executionTree/TernaryNode.h b/src/function/executionTree/TernaryNode.h
index be3b8e89b7e20f3dfc44f7de1a5edda2787bbf4f..3ede6cc6aa619234f80cbcba44f01f161bd3daf4 100644
--- a/src/function/executionTree/TernaryNode.h
+++ b/src/function/executionTree/TernaryNode.h
@@ -115,7 +115,7 @@ namespace gmshfem::function
       {
         bool isAllConstant = true;
         for(int i = 0; i < 3; ++i) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp critical
 #endif
           enclose._isConstant[i] = enclose._leaves.get(i)->isConstant(entity);
@@ -189,7 +189,7 @@ namespace gmshfem::function
             }
             return;
           });
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
         common::turn_into_tuple< 0, 3 >(
@@ -322,21 +322,21 @@ namespace gmshfem::function
     virtual void evaluate(OutputVector< typename T_Op::D_Scalar, T_Op::D_Degree > &values, const std::vector< scalar::Precision< typename T_Op::D_Scalar >, numa::allocator< scalar::Precision< typename T_Op::D_Scalar > > > &points, const std::vector< scalar::Precision< typename T_Op::D_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const override
     {
       if(values.size() < points.size() / 3) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
         FunctionAllocator< typename MathObject< typename T_Op::D_Scalar, T_Op::D_Degree >::Object >::initalizeNUMAMemory(points.size() / 3);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp master
 #endif
         values.resize(points.size() / 3);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
       }
       EfficientEvaluate evaluator;
       evaluator(*this, values, points, gaussPoints, elementType, entity);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 
 #pragma omp single
@@ -345,14 +345,14 @@ namespace gmshfem::function
         _tmpA.clear();
         _tmpA.shrink_to_fit();
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
         _tmpB.clear();
         _tmpB.shrink_to_fit();
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
diff --git a/src/function/executionTree/UnaryNode.h b/src/function/executionTree/UnaryNode.h
index a2f2c53fa5cfc61b4e9f2dfffa0919c5f2f13787..713944f14088410e989c3d518fbe6f74945b0b03 100644
--- a/src/function/executionTree/UnaryNode.h
+++ b/src/function/executionTree/UnaryNode.h
@@ -93,7 +93,7 @@ namespace gmshfem::function
         }
 
         enclose._leaves.get(0)->evaluate(*A, isConstant ? enclose._zeros : points, gaussPoints, elementType, entity);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
         if(isConstant) {
@@ -144,21 +144,21 @@ namespace gmshfem::function
     virtual void evaluate(OutputVector< typename T_Op::D_Scalar, T_Op::D_Degree > &values, const std::vector< scalar::Precision< typename T_Op::D_Scalar >, numa::allocator< scalar::Precision< typename T_Op::D_Scalar > > > &points, const std::vector< scalar::Precision< typename T_Op::D_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const override
     {
       if(values.size() < points.size() / 3) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
         FunctionAllocator< typename MathObject< typename T_Op::D_Scalar, T_Op::D_Degree >::Object >::initalizeNUMAMemory(points.size() / 3);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp master
 #endif
         values.resize(points.size() / 3);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
       }
       EfficientEvaluate evaluator;
       evaluator(*this, values, points, gaussPoints, elementType, entity);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 
 #pragma omp single
diff --git a/src/function/executionTree/binaryOperations.h b/src/function/executionTree/binaryOperations.h
index 52f32bb94499ce4271d66e0159e8853c65423dbf..74801538a3d4879c70a82ec0425187f741caa4c8 100644
--- a/src/function/executionTree/binaryOperations.h
+++ b/src/function/executionTree/binaryOperations.h
@@ -45,7 +45,7 @@ namespace gmshfem::function
     template< class T_A, class T_B >
     void operator()(OutputVector< T_Scalar, T_Degree > &values, const InputVector< T_Scalar, T_Degree, T_A > &a, const InputVector< T_Scalar, T_Degree, T_B > &b, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -84,7 +84,7 @@ namespace gmshfem::function
     template< class T_A, class T_B >
     void operator()(OutputVector< T_Scalar, T_Degree > &values, const InputVector< T_Scalar, T_Degree, T_A > &a, const InputVector< T_Scalar, T_Degree, T_B > &b, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -123,7 +123,7 @@ namespace gmshfem::function
     template< class T_A, class T_B >
     void operator()(OutputVector< T_Scalar, Degree::Degree1 > &values, const InputVector< T_Scalar, Degree::Degree1, T_A > &a, const InputVector< T_Scalar, Degree::Degree1, T_B > &b, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -158,7 +158,7 @@ namespace gmshfem::function
     template< class T_A, class T_B >
     void operator()(OutputVector< std::complex< T_ComplexScalar >, Degree::Degree1 > &values, const InputVector< std::complex< T_ComplexScalar >, Degree::Degree1, T_A > &a, const InputVector< std::complex< T_ComplexScalar >, Degree::Degree1, T_B > &b, const std::vector< scalar::Precision< std::complex< T_ComplexScalar > >, numa::allocator< scalar::Precision< std::complex< T_ComplexScalar > > > > &points, const std::vector< scalar::Precision< std::complex< T_ComplexScalar > > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -197,7 +197,7 @@ namespace gmshfem::function
     template< class T_A, class T_B >
     void operator()(OutputVector< T_Scalar, T_Degree > &values, const InputVector< T_Scalar, T_Degree, T_A > &a, const InputVector< T_Scalar, Degree::Degree0, T_B > &b, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -232,7 +232,7 @@ namespace gmshfem::function
     template< class T_A, class T_B >
     void operator()(OutputVector< T_Scalar, Degree::Degree4 > &values, const InputVector< T_Scalar, Degree::Degree4, T_A > &a, const InputVector< T_Scalar, Degree::Degree0, T_B > &b, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -278,7 +278,7 @@ namespace gmshfem::function
     template< class T_A, class T_B >
     void operator()(OutputVector< T_Scalar, Degree::Degree2 > &values, const InputVector< T_Scalar, Degree::Degree1, T_A > &a, const InputVector< T_Scalar, Degree::Degree1, T_B > &b, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -313,7 +313,7 @@ namespace gmshfem::function
     template< class T_A, class T_B >
     void operator()(OutputVector< T_Scalar, Degree::Degree4 > &values, const InputVector< T_Scalar, Degree::Degree2, T_A > &a, const InputVector< T_Scalar, Degree::Degree2, T_B > &b, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -360,7 +360,7 @@ namespace gmshfem::function
     template< class T_A, class T_B >
     void operator()(OutputVector< T_Scalar, T_Degree > &values, const InputVector< T_Scalar, T_Degree, T_A > &a, const InputVector< T_Scalar, T_Degree, T_B > &b, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -395,7 +395,7 @@ namespace gmshfem::function
     template< class T_A, class T_B >
     void operator()(OutputVector< T_Scalar, Degree::Degree0 > &values, const InputVector< T_Scalar, Degree::Degree0, T_A > &a, const InputVector< T_Scalar, Degree::Degree0, T_B > &b, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -430,7 +430,7 @@ namespace gmshfem::function
     template< class T_A, class T_B >
     void operator()(OutputVector< T_Scalar, Degree::Degree4 > &values, const InputVector< T_Scalar, Degree::Degree4, T_A > &a, const InputVector< T_Scalar, Degree::Degree4, T_B > &b, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -473,7 +473,7 @@ namespace gmshfem::function
     template< class T_A, class T_B >
     void operator()(OutputVector< T_Scalar, MathProduct< T_DegreeA, T_DegreeB >::ProductDegree > &values, const InputVector< T_Scalar, T_DegreeA, T_A > &a, const InputVector< T_Scalar, T_DegreeB, T_B > &b, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -508,7 +508,7 @@ namespace gmshfem::function
     template< class T_A, class T_B >
     void operator()(OutputVector< T_Scalar, MathProduct< Degree::Degree0, Degree::Degree4 >::ProductDegree > &values, const InputVector< T_Scalar, Degree::Degree0, T_A > &a, const InputVector< T_Scalar, Degree::Degree4, T_B > &b, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -547,7 +547,7 @@ namespace gmshfem::function
     template< class T_A, class T_B >
     void operator()(OutputVector< T_Scalar, MathProduct< Degree::Degree1, Degree::Degree1 >::ProductDegree > &values, const InputVector< T_Scalar, Degree::Degree1, T_A > &a, const InputVector< T_Scalar, Degree::Degree1, T_B > &b, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -577,7 +577,7 @@ namespace gmshfem::function
     template< class T_A, class T_B >
     void operator()(OutputVector< T_Scalar, MathProduct< Degree::Degree1, Degree::Degree2 >::ProductDegree > &values, const InputVector< T_Scalar, Degree::Degree1, T_A > &a, const InputVector< T_Scalar, Degree::Degree2, T_B > &b, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -612,7 +612,7 @@ namespace gmshfem::function
     template< class T_A, class T_B >
     void operator()(OutputVector< T_Scalar, MathProduct< Degree::Degree2, Degree::Degree2 >::ProductDegree > &values, const InputVector< T_Scalar, Degree::Degree2, T_A > &a, const InputVector< T_Scalar, Degree::Degree2, T_B > &b, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -647,7 +647,7 @@ namespace gmshfem::function
     template< class T_A, class T_B >
     void operator()(OutputVector< T_Scalar, MathProduct< Degree::Degree4, Degree::Degree0 >::ProductDegree > &values, const InputVector< T_Scalar, Degree::Degree4, T_A > &a, const InputVector< T_Scalar, Degree::Degree0, T_B > &b, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -688,7 +688,7 @@ namespace gmshfem::function
     template< class T_A, class T_B >
     void operator()(OutputVector< T_Scalar, MathProduct< Degree::Degree4, Degree::Degree2 >::ProductDegree > &values, const InputVector< T_Scalar, Degree::Degree4, T_A > &a, const InputVector< T_Scalar, Degree::Degree2, T_B > &b, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -729,7 +729,7 @@ namespace gmshfem::function
     template< class T_A, class T_B >
     void operator()(OutputVector< T_Scalar, MathProduct< Degree::Degree4, Degree::Degree4 >::ProductDegree > &values, const InputVector< T_Scalar, Degree::Degree4, T_A > &a, const InputVector< T_Scalar, Degree::Degree4, T_B > &b, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -785,7 +785,7 @@ namespace gmshfem::function
     template< class T_A, class T_B >
     void operator()(OutputVector< T_Scalar, T_Degree > &values, const InputVector< T_Scalar, T_Degree, T_A > &a, const InputVector< T_Scalar, T_Degree, T_B > &b, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
diff --git a/src/function/executionTree/fieldOperations.h b/src/function/executionTree/fieldOperations.h
index d6d9d99d1dc88d3be1ca8f8ccf1423703a6a6ee5..c57d6c8a7cf02a0b1b0247e5dee9848eb1d82df3 100644
--- a/src/function/executionTree/fieldOperations.h
+++ b/src/function/executionTree/fieldOperations.h
@@ -62,7 +62,7 @@ namespace gmshfem::function
     {
       // If the current model is not the same on which the field is defined, we use the point evaluation algorithm
       // otherwise, we can interpolate the field as usual.
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
@@ -74,18 +74,18 @@ namespace gmshfem::function
       }
       if(_pointEvaluation) {
         double verbose = 0.;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp master
 #endif
         {
           gmsh::option::getNumber("General.Verbosity", verbose);
           gmsh::option::setNumber("General.Verbosity", 0);
         }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         gmsh::model::mesh::rebuildElementCache(true);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < points.size() / 3; ++i) {
@@ -138,7 +138,7 @@ namespace gmshfem::function
           _fieldEvaluator(fieldExpression, &functions[0], &jacobians[0], &determinants[0], nbrDofsByElements);
           MathObject< T_Scalar, field::DegreeOfForm< field::NextForm< T_Form >::value >::value >::copy(values[i], fieldExpression * vecDof);
         }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp master
 #endif
         {
@@ -150,13 +150,13 @@ namespace gmshfem::function
         const unsigned int numThreads = omp::getNumThreads();
         const unsigned int myThreadID = omp::getThreadNum();
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         _nbrDofsByElements = field->getFunctionSpace()->getBasisFunctions(true, _functions, _fsIndex, gaussPoints, elementType, entity);
 
         std::vector< scalar::Precision< T_Scalar > > coordinates;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         field->getFunctionSpace()->getKeys(false, _typeKeys, _entityKeys, coordinates, elementType, entity);
@@ -165,7 +165,7 @@ namespace gmshfem::function
 
         const bool needOffset = !field->getFunctionSpace()->isConstantByElements();
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         _dofsValues.resize(_typeKeys.size());
@@ -173,39 +173,39 @@ namespace gmshfem::function
         field->getValues(_typeKeys, _entityKeys, _dofsValues, begin, end);
         const unsigned int nbrGaussPoints = gaussPoints.size() / 3;
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         scalar::copy(_gmshGaussPoints, gaussPoints);
 
         if(_fieldEvaluator.needJacobians()) {
           const domain::Domain domain = field->domain();
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
           gmsh::model::mesh::preallocateJacobians(elementType, nbrGaussPoints, true, true, false, _gmshJacobians, _gmshDeterminants, _gmshPoints, entity.second);
           gmsh::model::mesh::getJacobians(elementType, _gmshGaussPoints, _gmshJacobians, _gmshDeterminants, _gmshPoints, entity.second, myThreadID, numThreads);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 
 #pragma omp single
 #endif
           numa::copy(_jacobians, _gmshJacobians);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
           numa::copy(_determinants, _gmshDeterminants);
 
           if(domain.haveJacobiansModificators(entity)) {
             domain.applyJacobiansModificator(points, _determinants, _jacobians, entity);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
           }
         }
 
         Eigen::MatrixX< scalar::Precision< T_Scalar > > fieldExpression(_fieldEvaluator.size(), _nbrDofsByElements);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0U; i < nbrOfElements; ++i) {
@@ -217,7 +217,7 @@ namespace gmshfem::function
         }
       }
       // Restore the current model
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
@@ -273,7 +273,7 @@ namespace gmshfem::function
     {
       // If the current model is not the same on which the field is defined, we use the point evaluation algorithm
       // otherwise, we can interpolate the field as usual.
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
@@ -285,18 +285,18 @@ namespace gmshfem::function
       }
       if(_pointEvaluation) {
         double verbose = 0.;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp master
 #endif
         {
           gmsh::option::getNumber("General.Verbosity", verbose);
           gmsh::option::setNumber("General.Verbosity", 0);
         }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         gmsh::model::mesh::rebuildElementCache(true);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < points.size() / 3; ++i) {
@@ -351,7 +351,7 @@ namespace gmshfem::function
           Eigen::MatrixX< T_Scalar > tmp = fieldExpression * vecDof;
           MathObject< T_Scalar, field::DegreeOfCompoundForm< field::NextForm< T_Form >::value >::value >::copy(values[i], tmp);
         }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp master
 #endif
         {
@@ -363,13 +363,13 @@ namespace gmshfem::function
         const unsigned int numThreads = omp::getNumThreads();
         const unsigned int myThreadID = omp::getThreadNum();
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         _nbrDofsByElements = field->getFunctionSpace()->getBasisFunctions(true, _functions, _fsIndex, gaussPoints, elementType, entity);
 
         std::vector< scalar::Precision< T_Scalar > > coordinates;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         field->getFunctionSpace()->getKeys(false, _typeKeys, _entityKeys, coordinates, elementType, entity);
@@ -378,7 +378,7 @@ namespace gmshfem::function
 
         const bool needOffset = !field->getFunctionSpace()->isConstantByElements();
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         _dofsValues.resize(T_NumFields * _typeKeys.size());
@@ -386,32 +386,32 @@ namespace gmshfem::function
         field->getValues(_typeKeys, _entityKeys, _dofsValues, begin, end);
         const unsigned int nbrGaussPoints = gaussPoints.size() / 3;
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         scalar::copy(_gmshGaussPoints, gaussPoints);
 
         if(_fieldEvaluator.needJacobians()) {
           const domain::Domain domain = field->domain();
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
           gmsh::model::mesh::preallocateJacobians(elementType, nbrGaussPoints, true, true, false, _gmshJacobians, _gmshDeterminants, _gmshPoints, entity.second);
           gmsh::model::mesh::getJacobians(elementType, _gmshGaussPoints, _gmshJacobians, _gmshDeterminants, _gmshPoints, entity.second, myThreadID, numThreads);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 
 #pragma omp single
 #endif
           numa::copy(_jacobians, _gmshJacobians);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
           numa::copy(_determinants, _gmshDeterminants);
 
           if(domain.haveJacobiansModificators(entity)) {
             domain.applyJacobiansModificator(points, _determinants, _jacobians, entity);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
           }
@@ -419,7 +419,7 @@ namespace gmshfem::function
 
         Eigen::MatrixX< scalar::Precision< T_Scalar > > fieldExpression(_fieldEvaluator.size(), T_NumFields * _nbrDofsByElements);
         fieldExpression.setZero();
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0U; i < nbrOfElements; ++i) {
@@ -432,7 +432,7 @@ namespace gmshfem::function
         }
       }
       // Restore the current model
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
@@ -488,7 +488,7 @@ namespace gmshfem::function
     {
       // If the current model is not the same on which the field is defined, we use the point evaluation algorithm
       // otherwise, we can interpolate the field as usual.
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
@@ -500,18 +500,18 @@ namespace gmshfem::function
       }
       if(_pointEvaluation) {
         double verbose = 0.;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp master
 #endif
         {
           gmsh::option::getNumber("General.Verbosity", verbose);
           gmsh::option::setNumber("General.Verbosity", 0);
         }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         gmsh::model::mesh::rebuildElementCache(true);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < points.size() / 3; ++i) {
@@ -564,7 +564,7 @@ namespace gmshfem::function
           _fieldEvaluator(fieldExpression, &functions[0], &jacobians[0], &determinants[0], nbrDofsByElements);
           MathObject< T_Scalar, field::DegreeOfForm< T_Form >::value >::copy(values[i], fieldExpression * vecDof);
         }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp master
 #endif
         {
@@ -576,13 +576,13 @@ namespace gmshfem::function
         const unsigned int numThreads = omp::getNumThreads();
         const unsigned int myThreadID = omp::getThreadNum();
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         _nbrDofsByElements = field->getFunctionSpace()->getBasisFunctions(false, _functions, _fsIndex, gaussPoints, elementType, entity);
 
         std::vector< scalar::Precision< T_Scalar > > coordinates;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         field->getFunctionSpace()->getKeys(false, _typeKeys, _entityKeys, coordinates, elementType, entity);
@@ -591,7 +591,7 @@ namespace gmshfem::function
 
         const bool needOffset = !field->getFunctionSpace()->isConstantByElements();
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         _dofsValues.resize(_typeKeys.size());
@@ -599,39 +599,39 @@ namespace gmshfem::function
         field->getValues(_typeKeys, _entityKeys, _dofsValues, begin, end);
         const unsigned int nbrGaussPoints = gaussPoints.size() / 3;
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         scalar::copy(_gmshGaussPoints, gaussPoints);
 
         if(_fieldEvaluator.needJacobians()) {
           const domain::Domain domain = field->domain();
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
           gmsh::model::mesh::preallocateJacobians(elementType, nbrGaussPoints, true, true, false, _gmshJacobians, _gmshDeterminants, _gmshPoints, entity.second);
           gmsh::model::mesh::getJacobians(elementType, _gmshGaussPoints, _gmshJacobians, _gmshDeterminants, _gmshPoints, entity.second, myThreadID, numThreads);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 
 #pragma omp single
 #endif
           numa::copy(_jacobians, _gmshJacobians);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
           numa::copy(_determinants, _gmshDeterminants);
 
           if(domain.haveJacobiansModificators(entity)) {
             domain.applyJacobiansModificator(points, _determinants, _jacobians, entity);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
           }
         }
 
         Eigen::MatrixX< scalar::Precision< T_Scalar > > fieldExpression(_fieldEvaluator.size(), _nbrDofsByElements);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0U; i < nbrOfElements; ++i) {
@@ -643,7 +643,7 @@ namespace gmshfem::function
         }
       }
       // Restore the current model
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
@@ -700,7 +700,7 @@ namespace gmshfem::function
     {
       // If the current model is not the same on which the field is defined, we use the point evaluation algorithm
       // otherwise, we can interpolate the field as usual.
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
@@ -712,18 +712,18 @@ namespace gmshfem::function
       }
       if(_pointEvaluation) {
         double verbose = 0.;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp master
 #endif
         {
           gmsh::option::getNumber("General.Verbosity", verbose);
           gmsh::option::setNumber("General.Verbosity", 0);
         }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         gmsh::model::mesh::rebuildElementCache(true);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < points.size() / 3; ++i) {
@@ -778,7 +778,7 @@ namespace gmshfem::function
           Eigen::MatrixX< T_Scalar > tmp = fieldExpression * vecDof;
           MathObject< T_Scalar, field::DegreeOfCompoundForm< T_Form >::value >::copy(values[i], tmp);
         }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp master
 #endif
         {
@@ -790,13 +790,13 @@ namespace gmshfem::function
         const unsigned int numThreads = omp::getNumThreads();
         const unsigned int myThreadID = omp::getThreadNum();
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         _nbrDofsByElements = field->getFunctionSpace()->getBasisFunctions(false, _functions, _fsIndex, gaussPoints, elementType, entity);
 
         std::vector< scalar::Precision< T_Scalar > > coordinates;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         field->getFunctionSpace()->getKeys(false, _typeKeys, _entityKeys, coordinates, elementType, entity);
@@ -805,7 +805,7 @@ namespace gmshfem::function
 
         const bool needOffset = !field->getFunctionSpace()->isConstantByElements();
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         _dofsValues.resize(T_NumFields * _typeKeys.size());
@@ -813,32 +813,32 @@ namespace gmshfem::function
         field->getValues(_typeKeys, _entityKeys, _dofsValues, begin, end);
         const unsigned int nbrGaussPoints = gaussPoints.size() / 3;
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         scalar::copy(_gmshGaussPoints, gaussPoints);
 
         if(_fieldEvaluator.needJacobians()) {
           const domain::Domain domain = field->domain();
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
           gmsh::model::mesh::preallocateJacobians(elementType, nbrGaussPoints, true, true, false, _gmshJacobians, _gmshDeterminants, _gmshPoints, entity.second);
           gmsh::model::mesh::getJacobians(elementType, _gmshGaussPoints, _gmshJacobians, _gmshDeterminants, _gmshPoints, entity.second, myThreadID, numThreads);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 
 #pragma omp single
 #endif
           numa::copy(_jacobians, _gmshJacobians);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
           numa::copy(_determinants, _gmshDeterminants);
 
           if(domain.haveJacobiansModificators(entity)) {
             domain.applyJacobiansModificator(points, _determinants, _jacobians, entity);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
           }
@@ -846,7 +846,7 @@ namespace gmshfem::function
 
         Eigen::MatrixX< scalar::Precision< T_Scalar > > fieldExpression(_fieldEvaluator.size(), T_NumFields * _nbrDofsByElements);
         fieldExpression.setZero();
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0U; i < nbrOfElements; ++i) {
@@ -859,7 +859,7 @@ namespace gmshfem::function
         }
       }
       // Restore the current model
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
@@ -921,7 +921,7 @@ namespace gmshfem::function
     {
       // If the current model is not the same on which the field is defined, we use the point evaluation algorithm
       // otherwise, we can interpolate the field as usual.
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
@@ -933,18 +933,18 @@ namespace gmshfem::function
       }
       if(_pointEvaluation) {
         double verbose = 0.;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp master
 #endif
         {
           gmsh::option::getNumber("General.Verbosity", verbose);
           gmsh::option::setNumber("General.Verbosity", 0);
         }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         gmsh::model::mesh::rebuildElementCache(true);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < points.size() / 3; ++i) {
@@ -1001,7 +1001,7 @@ namespace gmshfem::function
           _fieldEvaluator(fieldExpression, &functions[0], &jacobians[0], &determinants[0], nbrDofsByElements);
           MathObject< scalar::ComplexPrecision< T_Scalar >, field::DegreeOfForm< T_Form >::value >::copy(values[i], fieldExpression * vecDof);
         }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp master
 #endif
         {
@@ -1013,13 +1013,13 @@ namespace gmshfem::function
         const unsigned int numThreads = omp::getNumThreads();
         const unsigned int myThreadID = omp::getThreadNum();
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         _nbrDofsByElements = field->getFunctionSpace()->getBasisFunctions(false, _functions, _fsIndex, gaussPoints, elementType, entity);
 
         std::vector< scalar::Precision< T_Scalar > > coordinates;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         field->getFunctionSpace()->getKeys(false, _typeKeys, _entityKeys, coordinates, elementType, entity);
@@ -1028,7 +1028,7 @@ namespace gmshfem::function
 
         const bool needOffset = !field->getFunctionSpace()->isConstantByElements();
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         _dofsValues.resize(_typeKeys.size());
@@ -1040,39 +1040,39 @@ namespace gmshfem::function
         module->getValues(_eigenTag, _typeKeys, _entityKeys, _dofsValues, begin, end);
         const unsigned int nbrGaussPoints = gaussPoints.size() / 3;
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         scalar::copy(_gmshGaussPoints, gaussPoints);
 
         if(_fieldEvaluator.needJacobians()) {
           const domain::Domain domain = field->domain();
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
           gmsh::model::mesh::preallocateJacobians(elementType, nbrGaussPoints, true, true, false, _gmshJacobians, _gmshDeterminants, _gmshPoints, entity.second);
           gmsh::model::mesh::getJacobians(elementType, _gmshGaussPoints, _gmshJacobians, _gmshDeterminants, _gmshPoints, entity.second, myThreadID, numThreads);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 
 #pragma omp single
 #endif
           numa::copy(_jacobians, _gmshJacobians);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
           numa::copy(_determinants, _gmshDeterminants);
 
           if(domain.haveJacobiansModificators(entity)) {
             domain.applyJacobiansModificator(points, _determinants, _jacobians, entity);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
           }
         }
 
         Eigen::MatrixX< scalar::Precision< T_Scalar > > fieldExpression(_fieldEvaluator.size(), _nbrDofsByElements);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0U; i < nbrOfElements; ++i) {
@@ -1084,7 +1084,7 @@ namespace gmshfem::function
         }
       }
       // Restore the current model
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
@@ -1141,7 +1141,7 @@ namespace gmshfem::function
     {
       // If the current model is not the same on which the field is defined, we use the point evaluation algorithm
       // otherwise, we can interpolate the field as usual.
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
@@ -1153,18 +1153,18 @@ namespace gmshfem::function
       }
       if(_pointEvaluation) {
         double verbose = 0.;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp master
 #endif
         {
           gmsh::option::getNumber("General.Verbosity", verbose);
           gmsh::option::setNumber("General.Verbosity", 0);
         }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         gmsh::model::mesh::rebuildElementCache(true);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < points.size() / 3; ++i) {
@@ -1223,7 +1223,7 @@ namespace gmshfem::function
           Eigen::MatrixX< scalar::ComplexPrecision< T_Scalar > > tmp = fieldExpression * vecDof;
           MathObject< scalar::ComplexPrecision< T_Scalar >, field::DegreeOfCompoundForm< T_Form >::value >::copy(values[i], tmp);
         }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp master
 #endif
         {
@@ -1235,13 +1235,13 @@ namespace gmshfem::function
         const unsigned int numThreads = omp::getNumThreads();
         const unsigned int myThreadID = omp::getThreadNum();
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         _nbrDofsByElements = field->getFunctionSpace()->getBasisFunctions(false, _functions, _fsIndex, gaussPoints, elementType, entity);
 
         std::vector< scalar::Precision< T_Scalar > > coordinates;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         field->getFunctionSpace()->getKeys(false, _typeKeys, _entityKeys, coordinates, elementType, entity);
@@ -1250,7 +1250,7 @@ namespace gmshfem::function
 
         const bool needOffset = !field->getFunctionSpace()->isConstantByElements();
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         _dofsValues.resize(T_NumFields * _typeKeys.size());
@@ -1262,32 +1262,32 @@ namespace gmshfem::function
         module->getValues(_eigenTag, _typeKeys, _entityKeys, _dofsValues, begin, end);
         const unsigned int nbrGaussPoints = gaussPoints.size() / 3;
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         scalar::copy(_gmshGaussPoints, gaussPoints);
 
         if(_fieldEvaluator.needJacobians()) {
           const domain::Domain domain = field->domain();
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
           gmsh::model::mesh::preallocateJacobians(elementType, nbrGaussPoints, true, true, false, _gmshJacobians, _gmshDeterminants, _gmshPoints, entity.second);
           gmsh::model::mesh::getJacobians(elementType, _gmshGaussPoints, _gmshJacobians, _gmshDeterminants, _gmshPoints, entity.second, myThreadID, numThreads);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 
 #pragma omp single
 #endif
           numa::copy(_jacobians, _gmshJacobians);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
           numa::copy(_determinants, _gmshDeterminants);
 
           if(domain.haveJacobiansModificators(entity)) {
             domain.applyJacobiansModificator(points, _determinants, _jacobians, entity);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 #endif
           }
@@ -1295,7 +1295,7 @@ namespace gmshfem::function
 
         Eigen::MatrixX< scalar::Precision< T_Scalar > > fieldExpression(_fieldEvaluator.size(), T_NumFields * _nbrDofsByElements);
         fieldExpression.setZero();
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0U; i < nbrOfElements; ++i) {
@@ -1308,7 +1308,7 @@ namespace gmshfem::function
         }
       }
       // Restore the current model
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
diff --git a/src/function/executionTree/novenaryOperations.h b/src/function/executionTree/novenaryOperations.h
index 15ee473e1482469af6842d871e6fdee5c3fe7970..b645c6d94000b06f83ae0e2d27a88f38c6da1b37 100644
--- a/src/function/executionTree/novenaryOperations.h
+++ b/src/function/executionTree/novenaryOperations.h
@@ -47,7 +47,7 @@ namespace gmshfem::function
                     const InputVector< T_Scalar, Degree::Degree0, T_I > &i,
                     const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto j = 0ULL; j < values.size(); ++j) {
@@ -90,7 +90,7 @@ namespace gmshfem::function
                     const InputVector< T_Scalar, Degree::Degree2, T_I > &i,
                     const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto j = 0ULL; j < values.size(); ++j) {
diff --git a/src/function/executionTree/nullaryOperations.h b/src/function/executionTree/nullaryOperations.h
index 950ebec5dfc6633424a1db5b918e6b72cab5ec4b..490b51ffa2d9a7ec0f94245d8f6ca7f0b6667b85 100644
--- a/src/function/executionTree/nullaryOperations.h
+++ b/src/function/executionTree/nullaryOperations.h
@@ -162,7 +162,7 @@ namespace gmshfem::function
       const std::vector< std::vector< typename MathObject< T_Scalar, T_Degree >::Object > > *pdata = (_dataSource ? _pdata : &_data);
       if(_method == 0) { // with _x and _y
         unsigned int i = 0, j = 0;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto n = 0ULL; n < values.size(); ++n) {
@@ -197,7 +197,7 @@ namespace gmshfem::function
             values[n] = scalar::Precision< T_Scalar >(1.) / ((x2 - x1) * (y2 - y1)) * (d11 * (x2 - xp) * (y2 - yp) + d12 * (x2 - xp) * (yp - y1) + d21 * (xp - x1) * (y2 - yp) + d22 * (xp - x1) * (yp - y1));
           }
           else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp critical
 #endif
             if(!warningOutside) {
@@ -213,7 +213,7 @@ namespace gmshfem::function
         const scalar::Precision< T_Scalar > coeff = 1. / (xStep * yStep);
         const scalar::Precision< T_Scalar > invXStep = 1. / xStep;
         const scalar::Precision< T_Scalar > invYStep = 1. / yStep;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto n = 0ULL; n < values.size(); ++n) {
@@ -247,7 +247,7 @@ namespace gmshfem::function
             values[n] = coeff * (d11 * (x2 - xp) * (y2 - yp) + d12 * (x2 - xp) * (yp - y1) + d21 * (xp - x1) * (y2 - yp) + d22 * (xp - x1) * (yp - y1));
           }
           else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp critical
 #endif
             if(!warningOutside) {
@@ -308,7 +308,7 @@ namespace gmshfem::function
 
     void operator()(OutputVector< T_Scalar, T_Degree > &values, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -354,7 +354,7 @@ namespace gmshfem::function
 
     void operator()(OutputVector< T_Scalar, Degree::Degree4 > &values, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -469,7 +469,7 @@ namespace gmshfem::function
       const std::vector< typename MathObject< T_Scalar, T_Degree >::Object > *pdata = (_dataSource ? _pdata : &_data);
       if(_method == 0) { // with _x and _y
         unsigned int i = 0;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto n = 0ULL; n < values.size(); ++n) {
@@ -492,7 +492,7 @@ namespace gmshfem::function
             values[n] = scalar::Precision< T_Scalar >(1.) / (x2 - x1) * (d1 * (x2 - xp) + d2 * (xp - x1));
           }
           else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp critical
 #endif
             if(!warningOutside) {
@@ -506,7 +506,7 @@ namespace gmshfem::function
         const scalar::Precision< T_Scalar > xStep = (_pMax - _pMin) / (_nx - 1);
         const scalar::Precision< T_Scalar > coeff = 1. / (xStep);
         const scalar::Precision< T_Scalar > invXStep = 1. / xStep;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto n = 0ULL; n < values.size(); ++n) {
@@ -526,7 +526,7 @@ namespace gmshfem::function
             values[n] = coeff * (d1 * (x2 - xp) + d2 * (xp - x1));
           }
           else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp critical
 #endif
             if(!warningOutside) {
@@ -596,7 +596,7 @@ namespace gmshfem::function
       }
 
       const unsigned int nbrOfElements = points.size() / gaussPoints.size();
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
@@ -607,7 +607,7 @@ namespace gmshfem::function
         scalar::move(_nodesCoord, gmshNodesCoord);
       }
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
@@ -617,7 +617,7 @@ namespace gmshfem::function
       const unsigned int nbrGaussPoints = gaussPoints.size() / 3;
 
       if(entity.first == 1) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0U; i < nbrOfElements; ++i) {
@@ -637,7 +637,7 @@ namespace gmshfem::function
         }
       }
       else if(entity.first == 2) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0U; i < nbrOfElements; ++i) {
@@ -694,7 +694,7 @@ namespace gmshfem::function
     void operator()(OutputVector< T_Scalar, Degree::Degree0 > &values, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
       // spherical angular coordinate (mathematics convention (r, theta, phi) with r is the radius, theta [0, 2pi] is the azimuthal angle is the xy plane and phi [0, pi] is the inclination polar angle)
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -763,7 +763,7 @@ namespace gmshfem::function
     {
       std::vector< double > viewval;
       double distance;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -806,7 +806,7 @@ namespace gmshfem::function
     {
       std::vector< double > viewval;
       double distance;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -849,7 +849,7 @@ namespace gmshfem::function
     {
       std::vector< double > viewval;
       double distance;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -892,7 +892,7 @@ namespace gmshfem::function
 
     void operator()(OutputVector< T_Scalar, Degree::Degree0 > &values, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -930,7 +930,7 @@ namespace gmshfem::function
 
     void operator()(OutputVector< T_Scalar, Degree::Degree0 > &values, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -983,7 +983,7 @@ namespace gmshfem::function
       }
 
       const unsigned int nbrOfElements = points.size() / gaussPoints.size();
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
@@ -994,7 +994,7 @@ namespace gmshfem::function
         scalar::move(_nodesCoord, gmshNodesCoord);
       }
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
@@ -1004,7 +1004,7 @@ namespace gmshfem::function
       const unsigned int nbrGaussPoints = gaussPoints.size() / 3;
 
       if(entity.first == 1) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0U; i < nbrOfElements; ++i) {
@@ -1024,18 +1024,18 @@ namespace gmshfem::function
         const unsigned int myThreadID = omp::getThreadNum();
         std::vector< double > gmshGaussPoints;
         scalar::copy(gmshGaussPoints, gaussPoints);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         gmsh::model::mesh::preallocateJacobians(elementType, nbrGaussPoints, true, true, false, _gmshJacobians, _gmshDeterminants, _gmshPoints, entity.second);
         gmsh::model::mesh::getJacobians(elementType, gmshGaussPoints, _gmshJacobians, _gmshDeterminants, _gmshPoints, entity.second, myThreadID, numThreads);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp barrier
 
 #pragma omp single
 #endif
         scalar::move(_jacobians, _gmshJacobians);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0U; i < nbrOfElements; ++i) {
@@ -1093,7 +1093,7 @@ namespace gmshfem::function
     void operator()(OutputVector< T_Scalar, Degree::Degree0 > &values, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
       // polar angular coordinate
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1287,7 +1287,7 @@ namespace gmshfem::function
       const std::vector< std::vector< std::vector< typename MathObject< T_Scalar, T_Degree >::Object > > > *pdata = (_dataSource ? _pdata : &_data);
       if(_method == 0) { // with _x _y and _z
         unsigned int i = 0, j = 0, k = 0;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto n = 0ULL; n < values.size(); ++n) {
@@ -1336,7 +1336,7 @@ namespace gmshfem::function
             values[n] = scalar::Precision< T_Scalar >(1.) / ((x2 - x1) * (y2 - y1) * (z2 - z1)) * (d111 * (x2 - xp) * (y2 - yp) * (z2 - zp) + d121 * (x2 - xp) * (yp - y1) * (z2 - zp) + d211 * (xp - x1) * (y2 - yp) * (z2 - zp) + d221 * (xp - x1) * (yp - y1) * (z2 - zp) + d112 * (x2 - xp) * (y2 - yp) * (zp - z1) + d122 * (x2 - xp) * (yp - y1) * (zp - z1) + d212 * (xp - x1) * (y2 - yp) * (zp - z1) + d222 * (xp - x1) * (yp - y1) * (zp - z1));
           }
           else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp critical
 #endif
             if(!warningOutside) {
@@ -1354,7 +1354,7 @@ namespace gmshfem::function
         const scalar::Precision< T_Scalar > invXStep = 1. / xStep;
         const scalar::Precision< T_Scalar > invYStep = 1. / yStep;
         const scalar::Precision< T_Scalar > invZStep = 1. / zStep;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto n = 0ULL; n < values.size(); ++n) {
@@ -1399,7 +1399,7 @@ namespace gmshfem::function
                                  d222 * (xp - x1) * (yp - y1) * (zp - z1));
           }
           else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp critical
 #endif
             if(!warningOutside) {
@@ -1474,7 +1474,7 @@ namespace gmshfem::function
 
     void operator()(OutputVector< T_Scalar, Degree::Degree0 > &values, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1512,7 +1512,7 @@ namespace gmshfem::function
 
     void operator()(OutputVector< T_Scalar, Degree::Degree0 > &values, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1550,7 +1550,7 @@ namespace gmshfem::function
 
     void operator()(OutputVector< T_Scalar, Degree::Degree0 > &values, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
diff --git a/src/function/executionTree/scalarTypeOperations.h b/src/function/executionTree/scalarTypeOperations.h
index 6dbcb22334951b969975e1706108cdcb20f4174f..e6e1402f38287b79408eb9d78e1ab9c027273c38 100644
--- a/src/function/executionTree/scalarTypeOperations.h
+++ b/src/function/executionTree/scalarTypeOperations.h
@@ -39,7 +39,7 @@ namespace gmshfem::function
 
     void operator()(OutputVector< scalar::Precision< T_Scalar >, T_Degree > &values, OutputVector< T_Scalar, T_Degree > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -68,7 +68,7 @@ namespace gmshfem::function
 
     void operator()(OutputVector< scalar::Precision< T_Scalar >, Degree::Degree1 > &values, OutputVector< T_Scalar, Degree::Degree1 > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -97,7 +97,7 @@ namespace gmshfem::function
 
     void operator()(OutputVector< scalar::Precision< T_Scalar >, Degree::Degree2 > &values, OutputVector< T_Scalar, Degree::Degree2 > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -126,7 +126,7 @@ namespace gmshfem::function
 
     void operator()(OutputVector< scalar::Precision< T_Scalar >, Degree::Degree4 > &values, OutputVector< T_Scalar, Degree::Degree4 > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -163,7 +163,7 @@ namespace gmshfem::function
 
     void operator()(OutputVector< scalar::Precision< T_Scalar >, T_Degree > &values, OutputVector< T_Scalar, T_Degree > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -192,7 +192,7 @@ namespace gmshfem::function
 
     void operator()(OutputVector< scalar::Precision< T_Scalar >, Degree::Degree1 > &values, OutputVector< T_Scalar, Degree::Degree1 > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -221,7 +221,7 @@ namespace gmshfem::function
 
     void operator()(OutputVector< scalar::Precision< T_Scalar >, Degree::Degree2 > &values, OutputVector< T_Scalar, Degree::Degree2 > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -250,7 +250,7 @@ namespace gmshfem::function
 
     void operator()(OutputVector< scalar::Precision< T_Scalar >, Degree::Degree4 > &values, OutputVector< T_Scalar, Degree::Degree4 > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -287,7 +287,7 @@ namespace gmshfem::function
 
     void operator()(OutputVector< scalar::ComplexPrecision< T_Scalar >, T_Degree > &values, OutputVector< T_Scalar, T_Degree > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -316,7 +316,7 @@ namespace gmshfem::function
 
     void operator()(OutputVector< scalar::ComplexPrecision< T_Scalar >, Degree::Degree1 > &values, OutputVector< T_Scalar, Degree::Degree1 > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -345,7 +345,7 @@ namespace gmshfem::function
 
     void operator()(OutputVector< scalar::ComplexPrecision< T_Scalar >, Degree::Degree2 > &values, OutputVector< T_Scalar, Degree::Degree2 > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -374,7 +374,7 @@ namespace gmshfem::function
 
     void operator()(OutputVector< scalar::ComplexPrecision< T_Scalar >, Degree::Degree4 > &values, OutputVector< T_Scalar, Degree::Degree4 > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
diff --git a/src/function/executionTree/ternaryOperations.h b/src/function/executionTree/ternaryOperations.h
index c92eede105d740808be6bf44ef732982abda3f49..1a3ee6a42366638c5a95d92e2576b85af5112d76 100644
--- a/src/function/executionTree/ternaryOperations.h
+++ b/src/function/executionTree/ternaryOperations.h
@@ -41,7 +41,7 @@ namespace gmshfem::function
                     const InputVector< T_Scalar, Degree::Degree0, T_C > &c,
                     const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto j = 0ULL; j < values.size(); ++j) {
@@ -79,7 +79,7 @@ namespace gmshfem::function
                     const InputVector< T_Scalar, Degree::Degree0, T_C > &c,
                     const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto j = 0ULL; j < values.size(); ++j) {
diff --git a/src/function/executionTree/unaryOperations.h b/src/function/executionTree/unaryOperations.h
index fd870c1d63c4e567c7259468d3d60385b2d19a22..44599fb5dfa00eaf50005a777c3ce99260e671e0 100644
--- a/src/function/executionTree/unaryOperations.h
+++ b/src/function/executionTree/unaryOperations.h
@@ -86,7 +86,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, T_Degree > &values, const InputVector< T_Scalar, T_Degree, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -125,7 +125,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, T_Degree > &values, const InputVector< T_Scalar, T_Degree, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -164,7 +164,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, T_Degree > &values, const InputVector< T_Scalar, T_Degree, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -203,7 +203,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, Degree::Degree0 > &values, const InputVector< T_Scalar, Degree::Degree1, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -233,7 +233,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< std::complex< T_ComplexScalar >, Degree::Degree0 > &values, const InputVector< std::complex< T_ComplexScalar >, Degree::Degree1, T_A > &a, const std::vector< scalar::Precision< std::complex< T_ComplexScalar > >, numa::allocator< scalar::Precision< std::complex< T_ComplexScalar > > > > &points, const std::vector< scalar::Precision< std::complex< T_ComplexScalar > > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -267,7 +267,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, T_Degree > &values, const InputVector< T_Scalar, T_Degree, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -306,7 +306,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, T_Degree > &values, const InputVector< T_Scalar, T_Degree, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -345,7 +345,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, T_Degree > &values, const InputVector< T_Scalar, T_Degree, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -384,7 +384,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, T_Degree > &values, const InputVector< T_Scalar, T_Degree, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -423,7 +423,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, T_Degree > &values, const InputVector< T_Scalar, T_Degree, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -458,7 +458,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< std::complex< T_ComplexScalar >, T_Degree > &values, const InputVector< std::complex< T_ComplexScalar >, T_Degree, T_A > &a, const std::vector< scalar::Precision< std::complex< T_ComplexScalar > >, numa::allocator< scalar::Precision< std::complex< T_ComplexScalar > > > > &points, const std::vector< scalar::Precision< std::complex< T_ComplexScalar > > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -503,7 +503,7 @@ namespace gmshfem::function
     void operator()(OutputVector< T_Scalar, T_Degree > &values, const InputVector< T_Scalar, T_Degree, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
       if(*_state) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < values.size(); ++i) {
@@ -511,7 +511,7 @@ namespace gmshfem::function
         }
       }
       else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < values.size(); ++i) {
@@ -559,7 +559,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, T_Degree > &values, const InputVector< T_Scalar, T_Degree, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -594,7 +594,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< std::complex< T_ComplexScalar >, Degree::Degree0 > &values, const InputVector< std::complex< T_ComplexScalar >, Degree::Degree0, T_A > &a, const std::vector< scalar::Precision< std::complex< T_ComplexScalar > >, numa::allocator< scalar::Precision< std::complex< T_ComplexScalar > > > > &points, const std::vector< scalar::Precision< std::complex< T_ComplexScalar > > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -629,7 +629,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< std::complex< T_ComplexScalar >, Degree::Degree1 > &values, const InputVector< std::complex< T_ComplexScalar >, Degree::Degree1, T_A > &a, const std::vector< scalar::Precision< std::complex< T_ComplexScalar > >, numa::allocator< scalar::Precision< std::complex< T_ComplexScalar > > > > &points, const std::vector< scalar::Precision< std::complex< T_ComplexScalar > > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -664,7 +664,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< std::complex< T_ComplexScalar >, Degree::Degree2 > &values, const InputVector< std::complex< T_ComplexScalar >, Degree::Degree2, T_A > &a, const std::vector< scalar::Precision< std::complex< T_ComplexScalar > >, numa::allocator< scalar::Precision< std::complex< T_ComplexScalar > > > > &points, const std::vector< scalar::Precision< std::complex< T_ComplexScalar > > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -703,7 +703,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, T_Degree > &values, const InputVector< T_Scalar, T_Degree, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -742,7 +742,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, T_Degree > &values, const InputVector< T_Scalar, T_Degree, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -787,13 +787,13 @@ namespace gmshfem::function
     {
 #ifdef HAVE_BESSEL
       double valr, vali;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
         int err = BesselJnComplex(_v, 1, double(std::real(a[i])), double(std::imag(a[i])), &valr, &vali);
         if(err != 0) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp critical
 #endif
           msg::error << "Issue with complex bessel function, error output: " << err << msg::endl;
@@ -846,13 +846,13 @@ namespace gmshfem::function
     {
 #ifdef HAVE_BESSEL
       double valr, vali;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
         int err = BesselJnComplex(_v, 1, double(std::real(a[i])), double(std::imag(a[i])), &valr, &vali);
         if(err != 0) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp critical
 #endif
           msg::error << "Issue with complex bessel function, error output: " << err << msg::endl;
@@ -909,13 +909,13 @@ namespace gmshfem::function
     {
 #ifdef HAVE_BESSEL
       double valr, vali;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
         int err = BesselYnComplex(_v, 1, double(std::real(a[i])), double(std::imag(a[i])), &valr, &vali);
         if(err != 0) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp critical
 #endif
           msg::error << "Issue with complex bessel function, error output: " << err << msg::endl;
@@ -968,13 +968,13 @@ namespace gmshfem::function
     {
 #ifdef HAVE_BESSEL
       double valr, vali;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
         int err = BesselYnComplex(_v, 1, double(std::real(a[i])), double(std::imag(a[i])), &valr, &vali);
         if(err != 0) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp critical
 #endif
           msg::error << "Issue with complex bessel function, error output:" << err << msg::endl;
@@ -1025,7 +1025,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, T_Degree > &values, const InputVector< T_Scalar, T_Degree, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1064,7 +1064,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, T_Degree > &values, const InputVector< T_Scalar, T_Degree, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1111,7 +1111,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, T_Degree > &values, const InputVector< T_Scalar, T_Degree, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1146,7 +1146,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< std::complex< T_ComplexScalar >, T_Degree > &values, const InputVector< std::complex< T_ComplexScalar >, T_Degree, T_A > &a, const std::vector< scalar::Precision< std::complex< T_ComplexScalar > >, numa::allocator< scalar::Precision< std::complex< T_ComplexScalar > > > > &points, const std::vector< scalar::Precision< std::complex< T_ComplexScalar > > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1185,7 +1185,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, T_Degree > &values, const InputVector< T_Scalar, T_Degree, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1220,7 +1220,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, Degree::Degree2 > &values, const InputVector< T_Scalar, Degree::Degree2, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1259,7 +1259,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, T_Degree > &values, const InputVector< T_Scalar, T_Degree, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1298,7 +1298,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, T_Degree > &values, const InputVector< T_Scalar, T_Degree, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1337,7 +1337,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, T_Degree > &values, const InputVector< T_Scalar, T_Degree, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1380,7 +1380,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, Degree::Degree0 > &values, const InputVector< T_Scalar, T_DegreeA, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1428,7 +1428,7 @@ namespace gmshfem::function
     void operator()(OutputVector< T_Scalar, Degree::Degree0 > &values, const InputVector< T_Scalar, Degree::Degree1, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
       if(_order == -1) { // Infinit norm
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1436,7 +1436,7 @@ namespace gmshfem::function
         }
       }
       else if(_order == 1) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1444,7 +1444,7 @@ namespace gmshfem::function
         }
       }
       else if(_order == 2) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1453,7 +1453,7 @@ namespace gmshfem::function
       }
       else {
         msg::warning << "Cannot apply the vectorial " + std::to_string(_order) + "-norm. The 2-norm is used" << msg::endl;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1497,7 +1497,7 @@ namespace gmshfem::function
     void operator()(OutputVector< T_Scalar, Degree::Degree0 > &values, const InputVector< T_Scalar, Degree::Degree2, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
       if(_order == -1) { // Infinit norm
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1505,7 +1505,7 @@ namespace gmshfem::function
         }
       }
       else if(_order == 1) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1513,7 +1513,7 @@ namespace gmshfem::function
         }
       }
       else if(_order == 2) { // Frobenius norm
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1522,7 +1522,7 @@ namespace gmshfem::function
       }
       else {
         msg::warning << "Cannot apply the tensorial " + std::to_string(_order) + "-norm. The 2-norm is used" << msg::endl;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1565,7 +1565,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, T_Degree > &values, const InputVector< T_Scalar, T_Degree, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1609,7 +1609,7 @@ namespace gmshfem::function
     void operator()(OutputVector< T_Scalar, T_Degree > &values, const InputVector< T_Scalar, T_Degree, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
       if(_pow == 1) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1617,7 +1617,7 @@ namespace gmshfem::function
         }
       }
       else if(_pow == 2) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1625,7 +1625,7 @@ namespace gmshfem::function
         }
       }
       else if(_pow == 3) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1633,7 +1633,7 @@ namespace gmshfem::function
         }
       }
       else if(_pow == 4) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1642,7 +1642,7 @@ namespace gmshfem::function
         }
       }
       else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1690,7 +1690,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, Degree::Degree0 > &values, const InputVector< T_Scalar, Degree::Degree1, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1720,7 +1720,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< std::complex< T_ComplexScalar >, Degree::Degree0 > &values, const InputVector< std::complex< T_ComplexScalar >, Degree::Degree1, T_A > &a, const std::vector< scalar::Precision< std::complex< T_ComplexScalar > >, numa::allocator< scalar::Precision< std::complex< T_ComplexScalar > > > > &points, const std::vector< scalar::Precision< std::complex< T_ComplexScalar > > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1754,7 +1754,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, T_Degree > &values, const InputVector< T_Scalar, T_Degree, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1789,7 +1789,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< std::complex< T_ComplexScalar >, T_Degree > &values, const InputVector< std::complex< T_ComplexScalar >, T_Degree, T_A > &a, const std::vector< scalar::Precision< std::complex< T_ComplexScalar > >, numa::allocator< scalar::Precision< std::complex< T_ComplexScalar > > > > &points, const std::vector< scalar::Precision< std::complex< T_ComplexScalar > > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1828,7 +1828,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, T_Degree > &values, const InputVector< T_Scalar, T_Degree, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1867,7 +1867,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, T_Degree > &values, const InputVector< T_Scalar, T_Degree, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1906,7 +1906,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, T_Degree > &values, const InputVector< T_Scalar, T_Degree, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1945,7 +1945,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, T_Degree > &values, const InputVector< T_Scalar, T_Degree, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -1984,7 +1984,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, T_Degree > &values, const InputVector< T_Scalar, T_Degree, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -2023,7 +2023,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, Degree::Degree0 > &values, const InputVector< T_Scalar, Degree::Degree2, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -2062,7 +2062,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, Degree::Degree0 > &values, const InputVector< T_Scalar, Degree::Degree1, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -2096,7 +2096,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, Degree::Degree0 > &values, const InputVector< T_Scalar, Degree::Degree2, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -2130,7 +2130,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, Degree::Degree0 > &values, const InputVector< T_Scalar, Degree::Degree2, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -2164,7 +2164,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, Degree::Degree0 > &values, const InputVector< T_Scalar, Degree::Degree2, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -2198,7 +2198,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, Degree::Degree0 > &values, const InputVector< T_Scalar, Degree::Degree1, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -2232,7 +2232,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, Degree::Degree0 > &values, const InputVector< T_Scalar, Degree::Degree2, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -2266,7 +2266,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, Degree::Degree0 > &values, const InputVector< T_Scalar, Degree::Degree2, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -2300,7 +2300,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, Degree::Degree0 > &values, const InputVector< T_Scalar, Degree::Degree2, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -2334,7 +2334,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, Degree::Degree0 > &values, const InputVector< T_Scalar, Degree::Degree1, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -2368,7 +2368,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, Degree::Degree0 > &values, const InputVector< T_Scalar, Degree::Degree2, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -2402,7 +2402,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, Degree::Degree0 > &values, const InputVector< T_Scalar, Degree::Degree2, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
@@ -2436,7 +2436,7 @@ namespace gmshfem::function
     template< class T_A >
     void operator()(OutputVector< T_Scalar, Degree::Degree0 > &values, const InputVector< T_Scalar, Degree::Degree2, T_A > &a, const std::vector< scalar::Precision< T_Scalar >, numa::allocator< scalar::Precision< T_Scalar > > > &points, const std::vector< scalar::Precision< T_Scalar > > &gaussPoints, const int elementType, const std::pair< int, int > &entity) const
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < values.size(); ++i) {
diff --git a/src/post/Integrate.cpp b/src/post/Integrate.cpp
index 56a9dc82d49f39820b20d1d592d7027a75047b62..1ed2d670f8d784e28fbef12b9cb3a737fef5c622 100644
--- a/src/post/Integrate.cpp
+++ b/src/post/Integrate.cpp
@@ -59,7 +59,7 @@ namespace gmshfem::post
 
         gmsh::model::mesh::preallocateJacobians(elementTypes[typeIndex], nbrOfGaussPoints, false, true, needPoints, gmshJacobians, gmshDeterminants, gmshPoints, itEntity->second);
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
         {
@@ -81,14 +81,14 @@ namespace gmshfem::post
           _function.evaluate(constant, 0., 0., 0., *itEntity);
         }
         std::vector< typename MathObject< T_Scalar, T_Degree >::Object, numa::allocator< typename MathObject< T_Scalar, T_Degree >::Object > > values; // Used only if _function isn't constant
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
         {
           const unsigned int myThreadID = omp::getThreadNum();
           if(_function.isConstant(*itEntity)) {
             const unsigned int nbrElements = determinants.size() / nbrOfGaussPoints;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
             for(auto j = 0U; j < nbrElements; ++j) {
@@ -100,7 +100,7 @@ namespace gmshfem::post
           else {
             const unsigned int nbrElements = determinants.size() / nbrOfGaussPoints;
             _function.evaluate(values, points, gaussPoints, elementTypes[typeIndex], *itEntity);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
             for(auto j = 0U; j < nbrElements; ++j) {
diff --git a/src/post/PostproMap.cpp b/src/post/PostproMap.cpp
index abef0549eb1a00269d708faabeb2a7e81cb0b669..7bb4b544ca387f1e21fad967b1bbbc360e1f2ece 100644
--- a/src/post/PostproMap.cpp
+++ b/src/post/PostproMap.cpp
@@ -34,7 +34,7 @@ namespace gmshfem::post
   {
     if(scalar::IsComplex< T_Scalar >::value) {
       // real part
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < gmshElementsTags.size(); ++i) {
@@ -42,13 +42,13 @@ namespace gmshfem::post
           gmshData[i].push_back(std::real(values[i * nbrNodesByElements + j]));
         }
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       gmsh::view::addModelData(tag, 2 * step, "", "ElementNodeData", gmshElementsTags, gmshData, time, 1, partition);
 
       // imaginary part
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < gmshElementsTags.size(); ++i) {
@@ -57,13 +57,13 @@ namespace gmshfem::post
           gmshData[i].push_back(std::imag(values[i * nbrNodesByElements + j]));
         }
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       gmsh::view::addModelData(tag, 2 * step + 1, "", "ElementNodeData", gmshElementsTags, gmshData, time, 1, partition);
     }
     else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < gmshElementsTags.size(); ++i) {
@@ -71,7 +71,7 @@ namespace gmshfem::post
           gmshData[i].push_back(std::real(values[i * nbrNodesByElements + j]));
         }
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       if(gmshData.size() != 0) {
@@ -85,7 +85,7 @@ namespace gmshfem::post
   {
     if(scalar::IsComplex< T_Scalar >::value) {
       // real part
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < gmshElementsTags.size(); ++i) {
@@ -96,13 +96,13 @@ namespace gmshfem::post
           }
         }
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       gmsh::view::addModelData(tag, 2 * step, "", "ElementNodeData", gmshElementsTags, gmshData, time, 3, partition);
 
       // imaginary part
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < gmshElementsTags.size(); ++i) {
@@ -114,13 +114,13 @@ namespace gmshfem::post
           }
         }
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       gmsh::view::addModelData(tag, 2 * step + 1, "", "ElementNodeData", gmshElementsTags, gmshData, time, 3, partition);
     }
     else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < gmshElementsTags.size(); ++i) {
@@ -131,7 +131,7 @@ namespace gmshfem::post
           }
         }
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       if(gmshData.size() != 0) {
@@ -145,7 +145,7 @@ namespace gmshfem::post
   {
     if(scalar::IsComplex< T_Scalar >::value) {
       // real part
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < gmshElementsTags.size(); ++i) {
@@ -158,13 +158,13 @@ namespace gmshfem::post
           }
         }
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       gmsh::view::addModelData(tag, 2 * step, "", "ElementNodeData", gmshElementsTags, gmshData, time, 9, partition);
 
       // imaginary part
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < gmshElementsTags.size(); ++i) {
@@ -178,13 +178,13 @@ namespace gmshfem::post
           }
         }
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       gmsh::view::addModelData(tag, 2 * step + 1, "", "ElementNodeData", gmshElementsTags, gmshData, time, 9, partition);
     }
     else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < gmshElementsTags.size(); ++i) {
@@ -197,7 +197,7 @@ namespace gmshfem::post
           }
         }
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       if(gmshData.size() != 0) {
@@ -241,7 +241,7 @@ namespace gmshfem::post
         gmsh::model::mesh::preallocateElementsByType(elementTypes[typeIndex], true, false, gmshElementsTags, gmshNodesTags, itEntity->second);
 
         std::vector< std::vector< double > > gmshData;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
         {
@@ -250,7 +250,7 @@ namespace gmshfem::post
 
           gmsh::model::mesh::getElementsByType(elementTypes[typeIndex], gmshElementsTags, gmshNodesTags, itEntity->second, myThreadID, numThreads);
           function.evaluate(values, coord, nodesCoord, elementTypes[typeIndex], *itEntity);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
           gmshData.resize(gmshElementsTags.size());
diff --git a/src/post/SaveField.cpp b/src/post/SaveField.cpp
index 0de59af2be3fa199e9d329b26fdfb7dda66ab8af..010f2c523f2e19ef947f2691544ab928769d76af 100644
--- a/src/post/SaveField.cpp
+++ b/src/post/SaveField.cpp
@@ -66,37 +66,37 @@ namespace gmshfem::post
     std::string dataType = "NodeData";
 
     if(scalar::IsComplex< T_Scalar >::value) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
       {
         // real part
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < data.size(); ++i) {
           gmshData[i] = {std::real(data[i])};
         }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         gmsh::view::addModelData(tag, 2 * _step, "", dataType, tags, gmshData, _time, 1, _partition);
 
         // imaginary part
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < data.size(); ++i) {
           gmshData[i] = {std::imag(data[i])};
         }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
         gmsh::view::addModelData(tag, 2 * _step + 1, "", dataType, tags, gmshData, _time, 1, _partition);
       }
     }
     else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
       for(auto i = 0ULL; i < data.size(); ++i) {
diff --git a/src/post/SaveFunction.cpp b/src/post/SaveFunction.cpp
index a9fcbbbe52a1cf5cd6e803d77a758c20a8718981..1a15f40f8d7469cb06a9c2842ceb83351ccc1ef1 100644
--- a/src/post/SaveFunction.cpp
+++ b/src/post/SaveFunction.cpp
@@ -31,7 +31,7 @@ namespace gmshfem::post
   {
     if(scalar::IsComplex< T_Scalar >::value) {
       // real part
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < gmshElementsTags.size(); ++i) {
@@ -39,13 +39,13 @@ namespace gmshfem::post
           gmshData[i].push_back(std::real(values[i * nbrNodesByElements + j]));
         }
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       gmsh::view::addModelData(tag, 2 * step, "", "ElementNodeData", gmshElementsTags, gmshData, time, 1, partition);
 
       // imaginary part
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < gmshElementsTags.size(); ++i) {
@@ -54,13 +54,13 @@ namespace gmshfem::post
           gmshData[i].push_back(std::imag(values[i * nbrNodesByElements + j]));
         }
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       gmsh::view::addModelData(tag, 2 * step + 1, "", "ElementNodeData", gmshElementsTags, gmshData, time, 1, partition);
     }
     else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < gmshElementsTags.size(); ++i) {
@@ -68,7 +68,7 @@ namespace gmshfem::post
           gmshData[i].push_back(std::real(values[i * nbrNodesByElements + j]));
         }
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       if(gmshData.size() != 0) {
@@ -82,7 +82,7 @@ namespace gmshfem::post
   {
     if(scalar::IsComplex< T_Scalar >::value) {
       // real part
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < gmshElementsTags.size(); ++i) {
@@ -93,13 +93,13 @@ namespace gmshfem::post
           }
         }
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       gmsh::view::addModelData(tag, 2 * step, "", "ElementNodeData", gmshElementsTags, gmshData, time, 3, partition);
 
       // imaginary part
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < gmshElementsTags.size(); ++i) {
@@ -111,13 +111,13 @@ namespace gmshfem::post
           }
         }
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       gmsh::view::addModelData(tag, 2 * step + 1, "", "ElementNodeData", gmshElementsTags, gmshData, time, 3, partition);
     }
     else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < gmshElementsTags.size(); ++i) {
@@ -128,7 +128,7 @@ namespace gmshfem::post
           }
         }
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       if(gmshData.size() != 0) {
@@ -142,7 +142,7 @@ namespace gmshfem::post
   {
     if(scalar::IsComplex< T_Scalar >::value) {
       // real part
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < gmshElementsTags.size(); ++i) {
@@ -155,13 +155,13 @@ namespace gmshfem::post
           }
         }
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       gmsh::view::addModelData(tag, 2 * step, "", "ElementNodeData", gmshElementsTags, gmshData, time, 9, partition);
 
       // imaginary part
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < gmshElementsTags.size(); ++i) {
@@ -175,13 +175,13 @@ namespace gmshfem::post
           }
         }
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       gmsh::view::addModelData(tag, 2 * step + 1, "", "ElementNodeData", gmshElementsTags, gmshData, time, 9, partition);
     }
     else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < gmshElementsTags.size(); ++i) {
@@ -194,7 +194,7 @@ namespace gmshfem::post
           }
         }
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       if(gmshData.size() != 0) {
@@ -239,7 +239,7 @@ namespace gmshfem::post
         gmsh::model::mesh::preallocateElementsByType(elementTypes[typeIndex], true, false, gmshElementsTags, gmshNodesTags, itEntity->second);
 
         std::vector< std::vector< double > > gmshData;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
         {
@@ -248,7 +248,7 @@ namespace gmshfem::post
 
           gmsh::model::mesh::getElementsByType(elementTypes[typeIndex], gmshElementsTags, gmshNodesTags, itEntity->second, myThreadID, numThreads);
           _function.evaluate(values, coord, nodesCoord, elementTypes[typeIndex], *itEntity);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
           gmshData.resize(gmshElementsTags.size());
diff --git a/src/problem/Formulation.cpp b/src/problem/Formulation.cpp
index 121dac852cd3b5923af0079987f6e37a60bc1e02..a9bc72f87063a073522534800bf0099bb10edfde 100644
--- a/src/problem/Formulation.cpp
+++ b/src/problem/Formulation.cpp
@@ -1502,7 +1502,7 @@ namespace gmshfem::problem
 
     gmsh::model::mesh::preallocateBarycenters(elementType, barycenters, tag);
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
     {
@@ -1534,12 +1534,12 @@ namespace gmshfem::problem
       std::vector< reorder::SortedEntity< 2 > > sortedEntities;
       sortedEntities.resize(nbrElements);
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
       {
         reorder::SortedEntity< 2 > se;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < nbrElements; ++i) {
@@ -1560,7 +1560,7 @@ namespace gmshfem::problem
       barycenters.clear();
       barycenters.shrink_to_fit();
       std::vector< std::size_t > degree(nbrElements);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
       for(auto i = 0ULL; i < nbrElements; ++i) {
@@ -1573,12 +1573,12 @@ namespace gmshfem::problem
       std::vector< reorder::SortedEntity< 3 > > sortedEntities;
       sortedEntities.resize(nbrElements);
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
       {
         reorder::SortedEntity< 3 > se;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < nbrElements; ++i) {
@@ -1602,7 +1602,7 @@ namespace gmshfem::problem
       barycenters.clear();
       barycenters.shrink_to_fit();
       std::vector< std::size_t > degree(nbrElements);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
       for(auto i = 0ULL; i < nbrElements; ++i) {
diff --git a/src/reorder/Hilbert.cpp b/src/reorder/Hilbert.cpp
index fe3291a81e6efb28c84b88c21d62c4a0729412c9..144095350325ebb22df2b9d93edcc1d6f91e5de4 100644
--- a/src/reorder/Hilbert.cpp
+++ b/src/reorder/Hilbert.cpp
@@ -193,7 +193,7 @@ namespace gmshfem::reorder
     uint64_t *histogramFull = new uint64_t[nbrBuckets * (numThreads + 1) + 1];
     uint64_t *histogram = &histogramFull[numThreads * nbrBuckets];
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
     {
@@ -201,14 +201,14 @@ namespace gmshfem::reorder
       uint64_t *histogramLocal = &histogramFull[myThreadID * nbrBuckets];
       memset(histogramLocal, 0, nbrBuckets * sizeof(uint64_t));
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
         memset(histogram, 0, (nbrBuckets + 1) * sizeof(uint64_t));
       }
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < nbrElements; ++i) {
@@ -225,7 +225,7 @@ namespace gmshfem::reorder
         histogramLocal[elements[i].h >> shift]++;
       }
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0U; i < nbrBuckets; ++i) {
@@ -237,7 +237,7 @@ namespace gmshfem::reorder
         }
       }
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
@@ -253,14 +253,14 @@ namespace gmshfem::reorder
         histogramLocal[i] += histogram[i];
       }
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < nbrElements; ++i) {
         elementsCopy[histogramLocal[elements[i].h >> shift]++] = elements[i];
       }
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for schedule(static, 1)
 #endif      
       for(auto i = 0U; i < nbrBuckets; ++i) {
diff --git a/src/reorder/RCM.cpp b/src/reorder/RCM.cpp
index a4ad4aa88960dfbc753b7b824f71c278313fba82..37db825de8e9bb0a0b8434607ca285dda35e14f8 100644
--- a/src/reorder/RCM.cpp
+++ b/src/reorder/RCM.cpp
@@ -32,13 +32,13 @@ namespace gmshfem::reorder
     std::vector< unsigned long long > Nmin;
     std::vector< unsigned long long > degree;
     std::vector< unsigned long long > newIndices;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
     {
       unsigned int numThreads = omp::getNumThreads();
       unsigned int myThreadID = omp::getThreadNum();
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
@@ -46,7 +46,7 @@ namespace gmshfem::reorder
         Nmin.resize(numThreads, 0);
         degree.resize(sorted.size(), 0);
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < sorted.size(); ++i) {
@@ -56,7 +56,7 @@ namespace gmshfem::reorder
           Nmin[myThreadID] = i;
         }
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
@@ -70,13 +70,13 @@ namespace gmshfem::reorder
       SortClass sortClass;
       sortClass._degree = degree.data();
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
         newIndices.resize(row[sorted.size()], 0);
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < sorted.size(); ++i) {
diff --git a/src/system/MatrixFactory.cpp b/src/system/MatrixFactory.cpp
index c5a098202af2c0ea17249fe55feda7efda937a5c..3e84ce608fb9f12c8cb526b483c72831de68c73d 100644
--- a/src/system/MatrixFactory.cpp
+++ b/src/system/MatrixFactory.cpp
@@ -29,17 +29,17 @@ namespace gmshfem::system
   {
     unsigned long long *newAJ = new unsigned long long[_ai[_size]];
     unsigned long long *newAI = new unsigned long long[_size + 1];
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < _size; ++i) {
         newAI[degree[i] - 1] = _ai[i + 1] - _ai[i];
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
@@ -51,7 +51,7 @@ namespace gmshfem::system
         }
         newAI[0] = 0;
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < _size; ++i) {
@@ -59,7 +59,7 @@ namespace gmshfem::system
           newAJ[newAI[degree[i] - 1] + j - _ai[i]] = degree[_aj[j]] - 1;
         }
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
       {
@@ -68,7 +68,7 @@ namespace gmshfem::system
         delete[] newAJ;
         delete[] newAI;
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
       for(auto i = 0ULL; i < _size; ++i) {
@@ -208,7 +208,7 @@ namespace gmshfem::system
 
     _aj = new unsigned long long[_ai[_size]];
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
     for(auto i = 0ULL; i < _size; ++i) {
@@ -285,18 +285,18 @@ namespace gmshfem::system
               HackComplex< scalar::Precision< T_Scalar > > *hack = reinterpret_cast< HackComplex< scalar::Precision< T_Scalar > > * >(&(*_module)[pos]);
               const scalar::Precision< T_Scalar > re = std::real(values[i * n + j]);
               const scalar::Precision< T_Scalar > im = std::imag(values[i * n + j]);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp atomic update
 #endif
               hack->r += re;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp atomic update
 #endif
               hack->i += im;
             }
             else {
               T_Scalar &update = (*_module)[pos];
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp atomic update
 #endif
               update += values[i * n + j];
@@ -308,11 +308,11 @@ namespace gmshfem::system
               const T_Scalar value = _valuesLC[indexJ[j].first] * values[i * n + j];
               const scalar::Precision< T_Scalar > re = std::real(value);
               const scalar::Precision< T_Scalar > im = std::imag(value);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp atomic update
 #endif
               hack->r += re;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp atomic update
 #endif
               hack->i += im;
@@ -320,7 +320,7 @@ namespace gmshfem::system
             else {
               T_Scalar &update = (*_module)[pos];
               const T_Scalar value = _valuesLC[indexJ[j].first] * values[i * n + j];
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp atomic update
 #endif
               update += value;
@@ -358,11 +358,11 @@ namespace gmshfem::system
               const T_Scalar value = std::conj(_valuesLC[indexI[i].first]) * values[i * n + j];
               const scalar::Precision< T_Scalar > re = std::real(value);
               const scalar::Precision< T_Scalar > im = std::imag(value);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp atomic update
 #endif
               hack->r += re;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp atomic update
 #endif
               hack->i += im;
@@ -370,7 +370,7 @@ namespace gmshfem::system
             else {
               T_Scalar &update = (*_module)[pos];
               const T_Scalar value = _valuesLC[indexI[i].first] * values[i * n + j];
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp atomic update
 #endif
               update += value;
@@ -382,11 +382,11 @@ namespace gmshfem::system
               const T_Scalar value = std::conj(_valuesLC[indexI[i].first]) * _valuesLC[indexJ[j].first] * values[i * n + j];
               const scalar::Precision< T_Scalar > re = std::real(value);
               const scalar::Precision< T_Scalar > im = std::imag(value);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp atomic update
 #endif
               hack->r += re;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp atomic update
 #endif
               hack->i += im;
@@ -394,7 +394,7 @@ namespace gmshfem::system
             else {
               T_Scalar &update = (*_module)[pos];
               const T_Scalar value = _valuesLC[indexI[i].first] * _valuesLC[indexJ[j].first] * values[i * n + j];
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp atomic update
 #endif
               update += value;
diff --git a/src/system/PetscInterface.h b/src/system/PetscInterface.h
index cae9b29333887e6a7da637529d3fcf2aadd1d297..b437afb93c630a61a435048b09e70739bda9312a 100644
--- a/src/system/PetscInterface.h
+++ b/src/system/PetscInterface.h
@@ -59,13 +59,13 @@ namespace gmshfem::system
       else {
         petscRow = (PetscInt *)std::malloc((size0 + 1) * sizeof(PetscInt));
         petscIndices = (PetscInt *)std::malloc(row[size0] * sizeof(PetscInt));
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
         for(auto i = 0ULL; i < size0 + 1; ++i) {
           petscRow[i] = row[i];
         }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
         for(auto i = 0ULL; i < row[size0]; ++i) {
@@ -91,7 +91,7 @@ namespace gmshfem::system
     // For solution
     void operator()(std::vector< T_Scalar1 > &values, T_Scalar1 *array, unsigned long long size)
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
       for(auto i = 0ULL; i < size; ++i) {
@@ -114,13 +114,13 @@ namespace gmshfem::system
       else {
         petscRow = (PetscInt *)std::malloc((size0 + 1) * sizeof(PetscInt));
         petscIndices = (PetscInt *)std::malloc(row[size0] * sizeof(PetscInt));
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
         for(auto i = 0ULL; i < size0 + 1; ++i) {
           petscRow[i] = row[i];
         }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
         for(auto i = 0ULL; i < row[size0]; ++i) {
@@ -147,7 +147,7 @@ namespace gmshfem::system
     // For solution
     void operator()(std::vector< T_Scalar1 > &values, T_Scalar1 *array, unsigned long long size)
     {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
       for(auto i = 0ULL; i < size; ++i) {
@@ -165,7 +165,7 @@ namespace gmshfem::system
         msg::warning << "Loss of efficiency can be noticed if you use a '" << scalar::Name< T_Scalar1 >::name << "' formulation with a '" << scalar::Name< T_Scalar2 >::name << "' PETSc" << msg::endl;
       }
       T_Scalar2 *petscValues = (T_Scalar2 *)std::malloc(row[size0] * sizeof(T_Scalar2));
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
       for(auto i = 0ULL; i < row[size0]; ++i) {
@@ -181,13 +181,13 @@ namespace gmshfem::system
       else {
         petscRow = (PetscInt *)std::malloc((size0 + 1) * sizeof(PetscInt));
         petscIndices = (PetscInt *)std::malloc(row[size0] * sizeof(PetscInt));
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
         for(auto i = 0ULL; i < size0 + 1; ++i) {
           petscRow[i] = row[i];
         }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
         for(auto i = 0ULL; i < row[size0]; ++i) {
@@ -213,7 +213,7 @@ namespace gmshfem::system
         msg::warning << "Loss of efficiency can be noticed if you use a '" << scalar::Name< T_Scalar1 >::name << "' formulation with a '" << scalar::Name< T_Scalar2 >::name << "' PETSc" << msg::endl;
       }
       T_Scalar2 *petscValues = (T_Scalar2 *)std::malloc(size * sizeof(T_Scalar2));
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
       for(auto i = 0ULL; i < size; ++i) {
@@ -229,7 +229,7 @@ namespace gmshfem::system
       if (sizeof(T_Scalar1) != sizeof(T_Scalar2)) {
         msg::warning << "The solution is cast because you use a '" << scalar::Name< T_Scalar1 >::name << "' formulation with a '" << scalar::Name< T_Scalar2 >::name << "' PETSc" << msg::endl;
       }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
       for(auto i = 0ULL; i < size; ++i) {
@@ -245,7 +245,7 @@ namespace gmshfem::system
     {
       msg::warning << "240 : Loss of efficiency can be noticed if you use a '" << scalar::Name< T_Scalar1 >::name << "' formulation with a '" << scalar::Name< T_Scalar2 >::name << "' PETSc" << msg::endl;
       T_Scalar2 *petscValues = (T_Scalar2 *)std::malloc(row[size0] * sizeof(T_Scalar2));
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
       for(auto i = 0ULL; i < row[size0]; ++i) {
@@ -261,13 +261,13 @@ namespace gmshfem::system
       else {
         petscRow = (PetscInt *)std::malloc((size0 + 1) * sizeof(PetscInt));
         petscIndices = (PetscInt *)std::malloc(row[size0] * sizeof(PetscInt));
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
         for(auto i = 0ULL; i < size0 + 1; ++i) {
           petscRow[i] = row[i];
         }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
         for(auto i = 0ULL; i < row[size0]; ++i) {
@@ -291,7 +291,7 @@ namespace gmshfem::system
     {
       msg::warning << "286 : Loss of efficiency can be noticed if you use a '" << scalar::Name< T_Scalar1 >::name << "' formulation with a '" << scalar::Name< T_Scalar2 >::name << "' PETSc" << msg::endl;
       T_Scalar2 *petscValues = (T_Scalar2 *)std::malloc(size * sizeof(T_Scalar2));
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
       for(auto i = 0ULL; i < size; ++i) {
@@ -305,7 +305,7 @@ namespace gmshfem::system
     void operator()(std::vector< T_Scalar1 > &values, T_Scalar2 *array, unsigned long long size)
     {
       msg::warning << "301 : The solution is cast because you use a '" << scalar::Name< T_Scalar1 >::name << "' formulation with a '" << scalar::Name< T_Scalar2 >::name << "' PETSc" << msg::endl;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
       for(auto i = 0ULL; i < size; ++i) {
@@ -321,7 +321,7 @@ namespace gmshfem::system
     {
       msg::warning << "316 : Loss of efficiency can be noticed if you use a '" << scalar::Name< T_Scalar1 >::name << "' formulation with a '" << scalar::Name< T_Scalar2 >::name << "' PETSc" << msg::endl;
       T_Scalar2 *petscValues = (T_Scalar2 *)std::malloc(row[size0] * sizeof(T_Scalar2));
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
       for(auto i = 0ULL; i < row[size0]; ++i) {
@@ -337,13 +337,13 @@ namespace gmshfem::system
       else {
         petscRow = (PetscInt *)std::malloc((size0 + 1) * sizeof(PetscInt));
         petscIndices = (PetscInt *)std::malloc(row[size0] * sizeof(PetscInt));
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
         for(auto i = 0ULL; i < size0 + 1; ++i) {
           petscRow[i] = row[i];
         }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
         for(auto i = 0ULL; i < row[size0]; ++i) {
@@ -367,7 +367,7 @@ namespace gmshfem::system
     {
       msg::warning << "362 : Loss of efficiency can be noticed if you use a '" << scalar::Name< T_Scalar1 >::name << "' formulation with a '" << scalar::Name< T_Scalar2 >::name << "' PETSc" << msg::endl;
       T_Scalar2 *petscValues = (T_Scalar2 *)std::malloc(size * sizeof(T_Scalar2));
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
       for(auto i = 0ULL; i < size; ++i) {
@@ -381,7 +381,7 @@ namespace gmshfem::system
     void operator()(std::vector< T_Scalar1 > &values, T_Scalar2 *array, unsigned long long size)
     {
       msg::warning << "377 : The solution is cast because you use a '" << scalar::Name< T_Scalar1 >::name << "' formulation with a '" << scalar::Name< T_Scalar2 >::name << "' PETSc" << msg::endl;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
       for(auto i = 0ULL; i < size; ++i) {
@@ -401,18 +401,18 @@ namespace gmshfem::system
       PetscInt *petscIndices = (PetscInt *)std::malloc((row[size0] * 4) * sizeof(PetscInt));
       T_Scalar2 *petscValues = (T_Scalar2 *)std::malloc((row[size0] * 4) * sizeof(T_Scalar2));
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
       {
         // row
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i <= size0; ++i) {
           petscRow[i] = 2 * row[i];
         }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 1ULL; i <= size0; ++i) {
@@ -420,7 +420,7 @@ namespace gmshfem::system
         }
 
         // indices
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < size0; ++i) {
@@ -431,7 +431,7 @@ namespace gmshfem::system
             petscIndices[petscRow[i] + j + row[i + 1] - 2 * row[i]] = indices[j] + size1;
           }
         }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < size0; ++i) {
@@ -444,7 +444,7 @@ namespace gmshfem::system
         }
 
         // values
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < size0; ++i) {
@@ -455,7 +455,7 @@ namespace gmshfem::system
             petscValues[petscRow[i] + j + row[i + 1] - 2 * row[i]] = (T_Scalar2(values[j].imag()) == 0. ? 0. : -T_Scalar2(values[j].imag()));
           }
         }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto i = 0ULL; i < size0; ++i) {
@@ -485,7 +485,7 @@ namespace gmshfem::system
       msg::warning << "479 : Loss of efficiency can be noticed if you use a '" << scalar::Name< T_Scalar1 >::name << "' formulation with a '" << scalar::Name< T_Scalar2 >::name << "' PETSc" << msg::endl;
 
       T_Scalar2 *petscValues = (T_Scalar2 *)std::malloc((size * 2) * sizeof(T_Scalar2));
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
       for(auto i = 0ULL; i < size; ++i) {
@@ -501,7 +501,7 @@ namespace gmshfem::system
     void operator()(std::vector< T_Scalar1 > &values, T_Scalar2 *array, unsigned long long size)
     {
       msg::warning << "497 : The solution is cast because you use a '" << scalar::Name< T_Scalar1 >::name << "' formulation with a '" << scalar::Name< T_Scalar2 >::name << "' PETSc" << msg::endl;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
       for(auto i = 0ULL; i < size / 2; ++i) {
diff --git a/src/system/VectorFactory.cpp b/src/system/VectorFactory.cpp
index 4a8781a8259949f427a7abc33a38800ccd1beec5..ce4f07d735ee7b187c05c0b80bbaf69b6e1d94ff 100644
--- a/src/system/VectorFactory.cpp
+++ b/src/system/VectorFactory.cpp
@@ -148,11 +148,11 @@ namespace gmshfem::system
         HackComplex< scalar::Precision< T_Scalar > > *hack = reinterpret_cast< HackComplex< scalar::Precision< T_Scalar > > * >(&_vec[index[i].first]);
         const scalar::Precision< T_Scalar > re = std::real(values[i]);
         const scalar::Precision< T_Scalar > im = std::imag(values[i]);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp atomic update
 #endif
         hack->r += re;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp atomic update
 #endif
         hack->i += im;
@@ -162,11 +162,11 @@ namespace gmshfem::system
         const T_Scalar value = std::conj(_valuesLC[index[i].first]) * values[i];
         const scalar::Precision< T_Scalar > re = std::real(value);
         const scalar::Precision< T_Scalar > im = std::imag(value);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp atomic update
 #endif
         hack->r += re;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp atomic update
 #endif
         hack->i += im;
@@ -179,13 +179,13 @@ namespace gmshfem::system
   {
     for(auto i = 0ULL; i < n; ++i) {
       if(index[i].second == dofs::AssembleType::Unknown || index[i].second == dofs::AssembleType::UnknownBubble || index[i].second == dofs::AssembleType::UnknownGlobal) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp atomic update
 #endif
         _vec[index[i].first] += values[i];
       }
       else if(index[i].second == dofs::AssembleType::Linked || index[i].second == dofs::AssembleType::LinkedBubble) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp atomic update
 #endif
         _vec[_indicesLC[index[i].first]] += _valuesLC[index[i].first] * values[i];
@@ -198,13 +198,13 @@ namespace gmshfem::system
   {
     for(auto i = 0ULL; i < n; ++i) {
       if(index[i].second == dofs::AssembleType::Unknown || index[i].second == dofs::AssembleType::UnknownBubble || index[i].second == dofs::AssembleType::UnknownGlobal) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp atomic update
 #endif
         _vec[index[i].first] += values[i];
       }
       else if(index[i].second == dofs::AssembleType::Linked || index[i].second == dofs::AssembleType::LinkedBubble) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp atomic update
 #endif
         _vec[_indicesLC[index[i].first]] += _valuesLC[index[i].first] * values[i];
@@ -223,11 +223,11 @@ namespace gmshfem::system
             HackComplex< scalar::Precision< T_Scalar > > *hack = reinterpret_cast< HackComplex< scalar::Precision< T_Scalar > > * >(&_vec[indexI[i].first]);
             const scalar::Precision< scalar::Precision< T_Scalar > > re = std::real(value);
             const scalar::Precision< scalar::Precision< T_Scalar > > im = std::imag(value);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp atomic update
 #endif
             hack->r += re;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp atomic update
 #endif
             hack->i += im;
@@ -241,11 +241,11 @@ namespace gmshfem::system
             HackComplex< scalar::Precision< T_Scalar > > *hack = reinterpret_cast< HackComplex< scalar::Precision< T_Scalar > > * >(&_vec[_indicesLC[indexI[i].first]]);
             const scalar::Precision< scalar::Precision< T_Scalar > > re = std::real(value);
             const scalar::Precision< scalar::Precision< T_Scalar > > im = std::imag(value);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp atomic update
 #endif
             hack->r += re;
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp atomic update
 #endif
             hack->i += im;
@@ -263,7 +263,7 @@ namespace gmshfem::system
         for(auto j = 0ULL; j < n; ++j) {
           if(indexJ[j].second == dofs::AssembleType::Fixed) {
             const double value = values[j * n + i] * _valuesDC[indexJ[j].first];
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp atomic update
 #endif
             _vec[indexI[i].first] += value;
@@ -274,7 +274,7 @@ namespace gmshfem::system
         for(auto j = 0ULL; j < n; ++j) {
           if(indexJ[j].second == dofs::AssembleType::Fixed) {
             const double value = values[j * n + i] * _valuesDC[indexJ[j].first];
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp atomic update
 #endif
             _vec[_indicesLC[indexI[i].first]] += _valuesLC[indexI[i].first] * value;
@@ -292,7 +292,7 @@ namespace gmshfem::system
         for(auto j = 0ULL; j < n; ++j) {
           if(indexJ[j].second == dofs::AssembleType::Fixed) {
             const float value = values[j * n + i] * _valuesDC[indexJ[j].first];
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp atomic update
 #endif
             _vec[indexI[i].first] += value;
@@ -303,7 +303,7 @@ namespace gmshfem::system
         for(auto j = 0ULL; j < n; ++j) {
           if(indexJ[j].second == dofs::AssembleType::Fixed) {
             const double value = values[j * n + i] * _valuesDC[indexJ[j].first];
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp atomic update
 #endif
             _vec[_indicesLC[indexI[i].first]] += _valuesLC[indexI[i].first] * value;
diff --git a/src/term/Assembler.cpp b/src/term/Assembler.cpp
index 5e44808e6b8d708fb131d5a18cca38d591a5e253..e1b484e76eb602e6b705f305daa0d6eb859c7350 100644
--- a/src/term/Assembler.cpp
+++ b/src/term/Assembler.cpp
@@ -118,7 +118,7 @@ namespace gmshfem::term
     for(const std::string &integrationType : integrationTypes) {
       unsigned int nbrOfGaussPoints = field::FunctionSpaceInterface< double >::GetGaussInfo(integrationType, elementType, integrationWeights, integrationPoints);
       gmsh::model::mesh::preallocateJacobians(elementType, nbrOfGaussPoints, needJacobians[integrationType], true, needPoints[integrationType], gmshJacobians, gmshDeterminants, gmshPoints, entity.second);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
       {
@@ -184,7 +184,7 @@ namespace gmshfem::term
           scalar::copy(gmshGaussPoints, gaussPoints);
 
           gmsh::model::mesh::preallocateJacobians(elementType, nbrOfGaussPoints, false, true, true, gmshJacobians, gmshDeterminants, gmshPoints, entity.second);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
           {
@@ -208,7 +208,7 @@ namespace gmshfem::term
             if(elements.jacobians(itIntegration->first) != nullptr) {
               jacobianJac->insert(jacobianJac->begin(), elements.jacobians(itIntegration->first)->begin(), elements.jacobians(itIntegration->first)->end());
             }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
             itIntegration->second[termIndex]->domain().applyJacobiansModificator(points, *determinantMod, *jacobianJac, entity);
@@ -230,7 +230,7 @@ namespace gmshfem::term
       }
       const std::string moduleName = A->getModule()->name();
       if(moduleName == "A") {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
         {
@@ -239,7 +239,7 @@ namespace gmshfem::term
           const unsigned long long nbrElements = determinants[terms[0]->tag()]->size() / terms[0]->nbrGaussPoints();
           Eigen::MatrixX< T_Scalar > A_e(nbrDofsByElementLhs, nbrDofsByElementRhs);
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
           for(auto elementIndex = 0ULL; elementIndex < nbrElements; ++elementIndex) {
@@ -262,7 +262,7 @@ namespace gmshfem::term
         if constexpr(scalar::IsComplex< T_Scalar >::value) {
           const scalar::ComplexPrecision< T_Scalar > im(0., 1.);
           const scalar::Precision< T_Scalar > frequency = static_cast< const system::AFrequencyModule< T_Scalar > * >(A->getModule())->getFrequency();
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
           {
@@ -272,7 +272,7 @@ namespace gmshfem::term
             Eigen::MatrixX< T_Scalar > A_e(nbrDofsByElementLhs, nbrDofsByElementRhs);
             Eigen::MatrixX< T_Scalar > A_e_tmp(nbrDofsByElementLhs, nbrDofsByElementRhs);
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
             for(auto elementIndex = 0ULL; elementIndex < nbrElements; ++elementIndex) {
@@ -306,12 +306,12 @@ namespace gmshfem::term
         }
       }
       else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
         for(auto idMatrix = 0ULL; idMatrix < moduleName.size(); ++idMatrix) {
           const char matrixType = moduleName[idMatrix];
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp single
 #endif
           A->getModule()->activate(matrixType);
@@ -320,7 +320,7 @@ namespace gmshfem::term
           const unsigned long long nbrElements = determinants[terms[0]->tag()]->size() / terms[0]->nbrGaussPoints();
           Eigen::MatrixX< T_Scalar > A_e(nbrDofsByElementLhs, nbrDofsByElementRhs);
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
           for(auto elementIndex = 0ULL; elementIndex < nbrElements; ++elementIndex) {
@@ -410,7 +410,7 @@ namespace gmshfem::term
           scalar::copy(gmshGaussPoints, gaussPoints);
 
           gmsh::model::mesh::preallocateJacobians(elementType, nbrOfGaussPoints, false, true, true, gmshJacobians, gmshDeterminants, gmshPoints, entity.second);
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
           {
@@ -434,7 +434,7 @@ namespace gmshfem::term
             if(elements.jacobians(integrationType) != nullptr) {
               jacobianJac->insert(jacobianJac->begin(), elements.jacobians(integrationType)->begin(), elements.jacobians(integrationType)->end());
             }
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
             termPtr->domain().applyJacobiansModificator(points, *determinantMod, *jacobianJac, entity);
@@ -451,7 +451,7 @@ namespace gmshfem::term
       }
 
       // assembly part
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
       {
@@ -459,7 +459,7 @@ namespace gmshfem::term
         const unsigned long long nbrElements = determinants[terms[0]->tag()]->size() / terms[0]->nbrGaussPoints();
         //Eigen::VectorX< T_Scalar > b_e(nbrDofsByElement);
         std::vector< Eigen::VectorX< T_Scalar > > b_es(b.size(), Eigen::VectorX< T_Scalar >(nbrDofsByElement));
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp for
 #endif
         for(auto elementIndex = 0ULL; elementIndex < nbrElements; ++elementIndex) {
diff --git a/src/term/BilinearTerm.cpp b/src/term/BilinearTerm.cpp
index 53ee2d1ea955b56d48004c8f368b64b1dc1c5bfa..1bafe7f3b62106d02eef92d9217af1362c58d6cb 100644
--- a/src/term/BilinearTerm.cpp
+++ b/src/term/BilinearTerm.cpp
@@ -179,7 +179,7 @@ namespace gmshfem::term
     }
     _fsIndexLhs = functionSpaces.functionSpaceOrientation(fs[0]->getGmshFemOrientationName());
     _fsLhs.resize(omp::getMaxThreads());
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
     {
@@ -200,7 +200,7 @@ namespace gmshfem::term
     _fsIndexRhs = functionSpaces.functionSpaceOrientation(fs[1]->getGmshFemOrientationName());
     _fsRhs.resize(omp::getMaxThreads());
     if(this->_field[0] != this->_field[1] || _isDerivativeLhs != _isDerivativeRhs) {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
       {
@@ -239,7 +239,7 @@ namespace gmshfem::term
       _fieldExpressionRhs->resize(omp::getMaxThreads(), Eigen::MatrixX< scalar::Precision< T_Scalar > >(_fieldEvaluatorRhs->size(), this->_field[1]->multiplicity() * this->_nbrOfDofsByElement[1]));
     }
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
     for(unsigned int i = 0; i < omp::getMaxThreads(); ++i) {
diff --git a/src/term/LinearTerm.cpp b/src/term/LinearTerm.cpp
index 62ea97a80c1fa7ebe1cf906bbaa1f852e63782ea..c35323b1ca06bd2c8c7f44f5fe261729a1215757 100644
--- a/src/term/LinearTerm.cpp
+++ b/src/term/LinearTerm.cpp
@@ -120,7 +120,7 @@ namespace gmshfem::term
       _lhsIsConstant = true;
     }
     else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
       {
@@ -138,7 +138,7 @@ namespace gmshfem::term
     // Instantiate tmp matrices
     // **************
     _fieldExpressionRhs.resize(omp::getMaxThreads(), Eigen::MatrixX< scalar::Precision< T_Scalar > >(_fieldEvaluatorRhs->size(), this->_fieldRhs->multiplicity() * this->_nbrOfDofsByElement));
-#ifdef _OPENMP  
+#ifdef HAVE_OPENMP  
   #pragma omp parallel for
 #endif
   for(unsigned int i = 0; i < omp::getMaxThreads(); ++i) {
diff --git a/src/term/Pattern.cpp b/src/term/Pattern.cpp
index 06c5e4c9257529ef77b82112d1e6eed213a84b1b..6d217a06589212c846082024fa4487ce853e951d 100644
--- a/src/term/Pattern.cpp
+++ b/src/term/Pattern.cpp
@@ -64,7 +64,7 @@ namespace gmshfem::term
       const std::vector< std::pair< unsigned long long, int > > *const indicesLhs = indices.indices(it->first->tag());
       const std::vector< std::pair< unsigned long long, int > > *const indicesRhs = indices.indices(it->second->tag());
 
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel for num_threads(omp::getMaxThreads())
 #endif
       for(auto elementIndex = 0ULL; elementIndex < indicesLhs->size() / nbrDofsByElement[0]; ++elementIndex) {
diff --git a/src/term/evaluator/EquationEvaluator.h b/src/term/evaluator/EquationEvaluator.h
index a9a48598b4ca2fe16f2ea7251a6164dca0b29e08..1e96a050f09bebf4b370d217e4d5489e647d95d7 100644
--- a/src/term/evaluator/EquationEvaluator.h
+++ b/src/term/evaluator/EquationEvaluator.h
@@ -160,7 +160,7 @@ namespace gmshfem::term::evaluator
         _left.evaluate(_values[0], 0., 0., 0., entity);
       }
       else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
         {
@@ -226,7 +226,7 @@ namespace gmshfem::term::evaluator
         _left.evaluate(_values[0], 0., 0., 0., entity);
       }
       else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
         {
@@ -292,7 +292,7 @@ namespace gmshfem::term::evaluator
         _left.evaluate(_values[0], 0., 0., 0., entity);
       }
       else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
         {
@@ -358,7 +358,7 @@ namespace gmshfem::term::evaluator
         _left.evaluate(_values[0], 0., 0., 0., entity);
       }
       else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
         {
@@ -424,7 +424,7 @@ namespace gmshfem::term::evaluator
         _left.evaluate(_values[0], 0., 0., 0., entity);
       }
       else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
         {
@@ -494,7 +494,7 @@ namespace gmshfem::term::evaluator
         _right.evaluate(_values[0], 0., 0., 0., entity);
       }
       else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
         {
@@ -560,7 +560,7 @@ namespace gmshfem::term::evaluator
         _right.evaluate(_values[0], 0., 0., 0., entity);
       }
       else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
         {
@@ -626,7 +626,7 @@ namespace gmshfem::term::evaluator
         _right.evaluate(_values[0], 0., 0., 0., entity);
       }
       else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
         {
@@ -692,7 +692,7 @@ namespace gmshfem::term::evaluator
         _right.evaluate(_values[0], 0., 0., 0., entity);
       }
       else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
         {
@@ -768,7 +768,7 @@ namespace gmshfem::term::evaluator
         _left.evaluate(_valuesLeft[0], 0., 0., 0., entity);
       }
       else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
         {
@@ -781,7 +781,7 @@ namespace gmshfem::term::evaluator
         _right.evaluate(_valuesRight[0], 0., 0., 0., entity);
       }
       else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
         {
@@ -880,7 +880,7 @@ namespace gmshfem::term::evaluator
         _left.evaluate(_values[0], 0., 0., 0., entity);
       }
       else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
         {
@@ -946,7 +946,7 @@ namespace gmshfem::term::evaluator
         _left.evaluate(_values[0], 0., 0., 0., entity);
       }
       else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
         {
@@ -1012,7 +1012,7 @@ namespace gmshfem::term::evaluator
         _left.evaluate(_values[0], 0., 0., 0., entity);
       }
       else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
         {
@@ -1078,7 +1078,7 @@ namespace gmshfem::term::evaluator
         _left.evaluate(_values[0], 0., 0., 0., entity);
       }
       else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
         {
@@ -1146,7 +1146,7 @@ namespace gmshfem::term::evaluator
         _left.evaluate(_values[0], 0., 0., 0., entity);
       }
       else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
         {
@@ -1214,7 +1214,7 @@ namespace gmshfem::term::evaluator
         _left.evaluate(_values[0], 0., 0., 0., entity);
       }
       else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
         {
@@ -1282,7 +1282,7 @@ namespace gmshfem::term::evaluator
         _left.evaluate(_values[0], 0., 0., 0., entity);
       }
       else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
         {
@@ -1360,7 +1360,7 @@ namespace gmshfem::term::evaluator
         _right.evaluate(_values[0], 0., 0., 0., entity);
       }
       else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
         {
@@ -1426,7 +1426,7 @@ namespace gmshfem::term::evaluator
         _right.evaluate(_values[0], 0., 0., 0., entity);
       }
       else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
         {
@@ -1492,7 +1492,7 @@ namespace gmshfem::term::evaluator
         _right.evaluate(_values[0], 0., 0., 0., entity);
       }
       else {
-#ifdef _OPENMP
+#ifdef HAVE_OPENMP
 #pragma omp parallel num_threads(omp::getMaxThreads())
 #endif
         {