diff --git a/tests/MeshDataBaseForTests.cpp b/tests/MeshDataBaseForTests.cpp
index 7afec8326a85a2b2e5f2e5f5ce46c2aae1f2b953..350dce0ffcd9562a3780b863be5d2872b191c7c7 100644
--- a/tests/MeshDataBaseForTests.cpp
+++ b/tests/MeshDataBaseForTests.cpp
@@ -22,6 +22,7 @@ MeshDataBaseForTests::MeshDataBaseForTests()
     CartesianMeshBuilder{TinyVector<3>{0, 1, 0}, TinyVector<3>{2, -1, 3}, TinyVector<3, size_t>{6, 7, 4}}.mesh());
 
   m_hybrid_2d_mesh = _buildHybrid2dMesh();
+  m_hybrid_3d_mesh = _buildHybrid3dMesh();
 }
 
 const MeshDataBaseForTests&
@@ -69,6 +70,12 @@ MeshDataBaseForTests::hybrid2DMesh() const
   return m_hybrid_2d_mesh;
 }
 
+std::shared_ptr<const Mesh<Connectivity<3>>>
+MeshDataBaseForTests::hybrid3DMesh() const
+{
+  return m_hybrid_3d_mesh;
+}
+
 std::shared_ptr<const Mesh<Connectivity<2>>>
 MeshDataBaseForTests::_buildHybrid2dMesh()
 {
@@ -240,3 +247,608 @@ $EndElements
   }
   return std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(GmshReader{filename}.mesh());
 }
+
+std::shared_ptr<const Mesh<Connectivity<3>>>
+MeshDataBaseForTests::_buildHybrid3dMesh()
+{
+  const std::string filename = std::filesystem::path{PUGS_BINARY_DIR}.append("tests").append("hybrid-3d.msh");
+  if (parallel::rank() == 0) {
+    std::ofstream fout(filename);
+    fout << R"($MeshFormat
+2.2 0 8
+$EndMeshFormat
+$PhysicalNames
+29
+0 40 "XMINYMINZMIN"
+0 41 "XMAXYMINZMIN"
+0 42 "XMINYMAXZMIN"
+0 43 "XMINYMAXZMAX"
+0 44 "XMINYMINZMAX"
+0 45 "XMAXYMINZMAX"
+0 47 "XMAXYMAXZMAX"
+0 51 "XMAXYMAXZMIN"
+1 28 "XMINZMIN"
+1 29 "XMINZMAX"
+1 30 "XMINYMAX"
+1 31 "XMINYMIN"
+1 32 "XMAXZMIN"
+1 33 "XMAXZMAX"
+1 34 "XMAXYMAX"
+1 35 "XMAXYMIN"
+1 36 "YMINZMIN"
+1 37 "YMINZMAX"
+1 38 "YMAXZMIN"
+1 39 "YMAXZMAX"
+2 22 "XMIN"
+2 23 "XMAX"
+2 24 "ZMAX"
+2 25 "ZMIN"
+2 26 "YMAX"
+2 27 "YMIN"
+3 52 "HEXAHEDRA"
+3 53 "PYRAMIDS"
+3 54 "TETRAHEDRA"
+$EndPhysicalNames
+$Nodes
+132
+1 0 0 0
+2 0.7 0 0
+3 0.8 1 0
+4 0 1 0
+5 1.3 0 0
+6 1.2 1 0
+7 0 1 1
+8 0 0 1
+9 0.7 0 1
+10 0.8 1 1
+11 1.3 0 1
+12 1.2 1 1
+13 2 0 0
+14 2 1 0
+15 2 0 1
+16 2 1 1
+17 0 0.7500000000003466 0
+18 0 0.5000000000020587 0
+19 0 0.2500000000010403 0
+20 0.7249999999999414 0.2499999999994139 0
+21 0.7499999999998691 0.4999999999986909 0
+22 0.7749999999999342 0.7499999999993421 0
+23 1.266666666666749 0.3333333333325114 0
+24 1.23333333333342 0.6666666666657952 0
+25 0.4000000000010956 1 0
+26 0.3499999999991014 0 0
+27 0.9999999999987851 0 0
+28 0 0.7500000000003466 1
+29 0 0.5000000000020587 1
+30 0 0.2500000000010403 1
+31 0.3499999999991014 0 1
+32 0.7249999999999414 0.2499999999994139 1
+33 0.7499999999998691 0.4999999999986909 1
+34 0.7749999999999342 0.7499999999993421 1
+35 0.4000000000010956 1 1
+36 0 1 0.3333333333333333
+37 0 1 0.6666666666666666
+38 0 0 0.3333333333333333
+39 0 0 0.6666666666666666
+40 0.7 0 0.3333333333333333
+41 0.7 0 0.6666666666666666
+42 0.8 1 0.3333333333333333
+43 0.8 1 0.6666666666666666
+44 0.9999999999987851 0 1
+45 1.266666666666749 0.3333333333325114 1
+46 1.23333333333342 0.6666666666657952 1
+47 1.3 0 0.3333333333333333
+48 1.3 0 0.6666666666666666
+49 1.2 1 0.3333333333333333
+50 1.2 1 0.6666666666666666
+51 1.630495225600928 0 0
+52 1.630495225600928 0 1
+53 2 0 0.499999999998694
+54 2 0.499999999998694 1
+55 2 0.499999999998694 0
+56 1.577708829260476 1 0
+57 1.577708829258421 1 1
+58 2 1 0.499999999998694
+59 0.3962888297748171 0.7916234456658734 0
+60 0.5588133104363274 0.5235006176096062 0
+61 0.48540966431768 0.3326974099157888 0
+62 0.2545886146233393 0.4410037988896774 0
+63 0.1888904744336107 0.2469120833355397 0
+64 0.9952483415423343 0.5591370310039079 0
+65 1.009047204638942 0.8056163713004886 0
+66 0.9496163434716166 0.3148023652713143 0
+67 0 0.7500000000003466 0.3333333333333333
+68 0 0.7500000000003466 0.6666666666666666
+69 0 0.5000000000020587 0.3333333333333333
+70 0 0.5000000000020587 0.6666666666666666
+71 0 0.2500000000010403 0.3333333333333333
+72 0 0.2500000000010403 0.6666666666666666
+73 0.3499999999991014 0 0.3333333333333333
+74 0.3499999999991014 0 0.6666666666666666
+75 0.7249999999999414 0.2499999999994139 0.3333333333333333
+76 0.7249999999999414 0.2499999999994139 0.6666666666666666
+77 0.7499999999998691 0.4999999999986909 0.3333333333333333
+78 0.7499999999998691 0.4999999999986909 0.6666666666666666
+79 0.7749999999999342 0.7499999999993421 0.3333333333333333
+80 0.7749999999999342 0.7499999999993421 0.6666666666666666
+81 0.4000000000010956 1 0.3333333333333333
+82 0.4000000000010956 1 0.6666666666666666
+83 0.3962888297748171 0.7916234456658734 1
+84 0.5588133104363274 0.5235006176096062 1
+85 0.48540966431768 0.3326974099157888 1
+86 0.2545886146233393 0.4410037988896774 1
+87 0.1888904744336107 0.2469120833355397 1
+88 0.9999999999987851 0 0.3333333333333333
+89 0.9999999999987851 0 0.6666666666666666
+90 1.266666666666749 0.3333333333325114 0.3333333333333333
+91 1.266666666666749 0.3333333333325114 0.6666666666666666
+92 1.23333333333342 0.6666666666657952 0.3333333333333333
+93 1.23333333333342 0.6666666666657952 0.6666666666666666
+94 0.9952483415423343 0.5591370310039079 1
+95 1.009047204638942 0.8056163713004886 1
+96 0.9496163434716166 0.3148023652713143 1
+97 1.694794622046484 0 0.4999999999999999
+98 2 0.667671095007814 0.3323289049917093
+99 2 0.3400507564846686 0.5028063863717376
+100 2 0.7126842003874787 0.6873157996117608
+101 1.638498486684473 0.2509397074193687 1
+102 1.619828855633658 0.6299359336032246 1
+103 1.638498486684556 0.2509397074193682 0
+104 1.619828855634071 0.629935933603235 0
+105 1.657210007479593 1 0.5000000082593427
+106 0.3962888297748171 0.7916234456658734 0.3333333333333333
+107 0.3962888297748171 0.7916234456658734 0.6666666666666666
+108 0.5588133104363274 0.5235006176096062 0.3333333333333333
+109 0.5588133104363274 0.5235006176096062 0.6666666666666666
+110 0.48540966431768 0.3326974099157888 0.3333333333333333
+111 0.48540966431768 0.3326974099157888 0.6666666666666666
+112 0.2545886146233393 0.4410037988896774 0.3333333333333333
+113 0.2545886146233393 0.4410037988896774 0.6666666666666666
+114 0.1888904744336107 0.2469120833355397 0.3333333333333333
+115 0.1888904744336107 0.2469120833355397 0.6666666666666666
+116 0.9952483415423343 0.5591370310039079 0.3333333333333333
+117 0.9952483415423343 0.5591370310039079 0.6666666666666666
+118 1.009047204638942 0.8056163713004886 0.3333333333333333
+119 1.009047204638942 0.8056163713004886 0.6666666666666666
+120 0.9496163434716166 0.3148023652713143 0.3333333333333333
+121 0.9496163434716166 0.3148023652713143 0.6666666666666666
+122 1.605048220195734 0.3680048220187183 0.5
+123 1.521560099439846 0.6946560099431293 0.4999999999999999
+124 1.449585918125941 0.1832919251455124 0.1666666666666667
+125 1.449585918125941 0.1832919251455124 0.4999999999999999
+126 1.449585918125941 0.1832919251455124 0.8333333333333333
+127 1.416252584792843 0.5166252584784292 0.1666666666666667
+128 1.416252584792843 0.5166252584784292 0.4999999999999999
+129 1.416252584792843 0.5166252584784292 0.8333333333333333
+130 1.382919251459699 0.8499585918121965 0.1666666666666667
+131 1.382919251459699 0.8499585918121965 0.4999999999999999
+132 1.382919251459699 0.8499585918121965 0.8333333333333333
+$EndNodes
+$Elements
+384
+1 15 2 40 1 1
+2 15 2 42 4 4
+3 15 2 43 7 7
+4 15 2 44 8 8
+5 15 2 41 27 13
+6 15 2 51 28 14
+7 15 2 45 29 15
+8 15 2 47 30 16
+9 1 2 28 1 4 17
+10 1 2 28 1 17 18
+11 1 2 28 1 18 19
+12 1 2 28 1 19 1
+13 1 2 38 4 3 25
+14 1 2 38 4 25 4
+15 1 2 38 5 6 3
+16 1 2 36 6 1 26
+17 1 2 36 6 26 2
+18 1 2 36 7 2 27
+19 1 2 36 7 27 5
+20 1 2 29 12 7 28
+21 1 2 29 12 28 29
+22 1 2 29 12 29 30
+23 1 2 29 12 30 8
+24 1 2 37 13 8 31
+25 1 2 37 13 31 9
+26 1 2 39 15 10 35
+27 1 2 39 15 35 7
+28 1 2 30 17 4 36
+29 1 2 30 17 36 37
+30 1 2 30 17 37 7
+31 1 2 31 18 1 38
+32 1 2 31 18 38 39
+33 1 2 31 18 39 8
+34 1 2 37 35 9 44
+35 1 2 37 35 44 11
+36 1 2 39 37 12 10
+37 1 2 36 49 5 51
+38 1 2 36 49 51 13
+39 1 2 37 50 11 52
+40 1 2 37 50 52 15
+41 1 2 35 51 13 53
+42 1 2 35 51 53 15
+43 1 2 33 52 15 54
+44 1 2 33 52 54 16
+45 1 2 32 53 13 55
+46 1 2 32 53 55 14
+47 1 2 38 54 14 56
+48 1 2 38 54 56 6
+49 1 2 39 55 12 57
+50 1 2 39 55 57 16
+51 1 2 34 56 14 58
+52 1 2 34 56 58 16
+53 2 2 25 2 20 2 27
+54 2 2 25 2 23 27 5
+55 2 2 25 2 64 23 24
+56 2 2 25 2 22 65 3
+57 2 2 25 2 23 64 66
+58 2 2 25 2 65 6 3
+59 2 2 25 2 64 21 66
+60 2 2 25 2 65 24 6
+61 2 2 25 2 24 65 64
+62 2 2 25 2 65 22 64
+63 2 2 25 2 66 20 27
+64 2 2 25 2 21 64 22
+65 2 2 25 2 20 66 21
+66 2 2 25 2 23 66 27
+67 2 2 24 54 32 9 44
+68 2 2 24 54 45 44 11
+69 2 2 24 54 94 45 46
+70 2 2 24 54 34 95 10
+71 2 2 24 54 45 94 96
+72 2 2 24 54 95 12 10
+73 2 2 24 54 94 33 96
+74 2 2 24 54 95 46 12
+75 2 2 24 54 46 95 94
+76 2 2 24 54 95 34 94
+77 2 2 24 54 96 32 44
+78 2 2 24 54 33 94 34
+79 2 2 24 54 32 96 33
+80 2 2 24 54 45 96 44
+81 2 2 27 55 51 5 47
+82 2 2 27 55 48 11 52
+83 2 2 27 55 51 97 13
+84 2 2 27 55 97 53 13
+85 2 2 27 55 97 52 15
+86 2 2 27 55 15 53 97
+87 2 2 27 55 97 47 48
+88 2 2 27 55 47 97 51
+89 2 2 27 55 48 52 97
+90 2 2 23 56 53 99 13
+91 2 2 23 56 99 55 13
+92 2 2 23 56 55 98 14
+93 2 2 23 56 98 58 14
+94 2 2 23 56 99 53 15
+95 2 2 23 56 54 99 15
+96 2 2 23 56 100 54 16
+97 2 2 23 56 58 100 16
+98 2 2 23 56 100 99 54
+99 2 2 23 56 55 99 98
+100 2 2 23 56 100 58 98
+101 2 2 23 56 98 99 100
+102 2 2 24 57 45 101 11
+103 2 2 24 57 101 52 11
+104 2 2 24 57 46 12 57
+105 2 2 24 57 52 101 15
+106 2 2 24 57 101 54 15
+107 2 2 24 57 54 102 16
+108 2 2 24 57 102 57 16
+109 2 2 24 57 46 102 45
+110 2 2 24 57 101 45 102
+111 2 2 24 57 102 46 57
+112 2 2 24 57 102 54 101
+113 2 2 25 58 23 103 5
+114 2 2 25 58 103 51 5
+115 2 2 25 58 24 6 56
+116 2 2 25 58 51 103 13
+117 2 2 25 58 103 55 13
+118 2 2 25 58 55 104 14
+119 2 2 25 58 104 56 14
+120 2 2 25 58 24 104 23
+121 2 2 25 58 103 23 104
+122 2 2 25 58 104 24 56
+123 2 2 25 58 104 55 103
+124 2 2 26 59 6 49 56
+125 2 2 26 59 50 12 57
+126 2 2 26 59 56 105 14
+127 2 2 26 59 105 58 14
+128 2 2 26 59 105 57 16
+129 2 2 26 59 16 58 105
+130 2 2 26 59 105 49 50
+131 2 2 26 59 49 105 56
+132 2 2 26 59 50 57 105
+133 3 2 25 1 59 60 21 22
+134 3 2 25 1 22 3 25 59
+135 3 2 25 1 18 62 59 17
+136 3 2 25 1 25 4 17 59
+137 3 2 25 1 62 63 26 61
+138 3 2 25 1 60 61 20 21
+139 3 2 25 1 62 18 19 63
+140 3 2 25 1 26 2 20 61
+141 3 2 25 1 61 60 59 62
+142 3 2 25 1 19 1 26 63
+143 3 2 22 19 4 17 67 36
+144 3 2 22 19 36 67 68 37
+145 3 2 22 19 37 68 28 7
+146 3 2 22 19 17 18 69 67
+147 3 2 22 19 67 69 70 68
+148 3 2 22 19 68 70 29 28
+149 3 2 22 19 18 19 71 69
+150 3 2 22 19 69 71 72 70
+151 3 2 22 19 70 72 30 29
+152 3 2 22 19 19 1 38 71
+153 3 2 22 19 71 38 39 72
+154 3 2 22 19 72 39 8 30
+155 3 2 27 23 1 26 73 38
+156 3 2 27 23 38 73 74 39
+157 3 2 27 23 39 74 31 8
+158 3 2 27 23 26 2 40 73
+159 3 2 27 23 73 40 41 74
+160 3 2 27 23 74 41 9 31
+161 3 2 26 31 3 25 81 42
+162 3 2 26 31 42 81 82 43
+163 3 2 26 31 43 82 35 10
+164 3 2 26 31 25 4 36 81
+165 3 2 26 31 81 36 37 82
+166 3 2 26 31 82 37 7 35
+167 3 2 24 32 83 84 33 34
+168 3 2 24 32 34 10 35 83
+169 3 2 24 32 29 86 83 28
+170 3 2 24 32 35 7 28 83
+171 3 2 24 32 86 87 31 85
+172 3 2 24 32 84 85 32 33
+173 3 2 24 32 86 29 30 87
+174 3 2 24 32 31 9 32 85
+175 3 2 24 32 85 84 83 86
+176 3 2 24 32 30 8 31 87
+177 3 2 27 45 2 27 88 40
+178 3 2 27 45 40 88 89 41
+179 3 2 27 45 41 89 44 9
+180 3 2 27 45 27 5 47 88
+181 3 2 27 45 88 47 48 89
+182 3 2 27 45 89 48 11 44
+183 3 2 26 53 6 3 42 49
+184 3 2 26 53 49 42 43 50
+185 3 2 26 53 50 43 10 12
+186 4 2 54 3 90 103 127 122
+187 4 2 54 3 101 91 129 122
+188 4 2 54 3 103 90 124 122
+189 4 2 54 3 126 91 101 122
+190 4 2 54 3 92 104 130 123
+191 4 2 54 3 102 93 132 123
+192 4 2 54 3 129 45 91 101
+193 4 2 54 3 90 23 127 103
+194 4 2 54 3 45 126 91 101
+195 4 2 54 3 23 90 124 103
+196 4 2 54 3 92 127 104 123
+197 4 2 54 3 93 102 129 123
+198 4 2 54 3 125 90 91 122
+199 4 2 54 3 91 90 128 122
+200 4 2 54 3 104 127 103 122
+201 4 2 54 3 101 129 102 122
+202 4 2 54 3 127 122 104 123
+203 4 2 54 3 102 122 129 123
+204 4 2 54 3 103 124 97 122
+205 4 2 54 3 126 101 97 122
+206 4 2 54 3 132 57 102 123
+207 4 2 54 3 130 104 56 123
+208 4 2 54 3 46 93 132 102
+209 4 2 54 3 92 24 130 104
+210 4 2 54 3 129 93 46 102
+211 4 2 54 3 92 127 24 104
+212 4 2 54 3 102 100 54 122
+213 4 2 54 3 100 122 102 123
+214 4 2 54 3 54 101 102 122
+215 4 2 54 3 104 103 55 122
+216 4 2 54 3 51 103 124 97
+217 4 2 54 3 126 101 52 97
+218 4 2 54 3 105 130 56 123
+219 4 2 54 3 132 105 57 123
+220 4 2 54 3 98 104 55 122
+221 4 2 54 3 93 92 131 123
+222 4 2 54 3 93 128 92 123
+223 4 2 54 3 51 124 47 97
+224 4 2 54 3 52 48 126 97
+225 4 2 54 3 54 100 99 122
+226 4 2 54 3 104 122 98 123
+227 4 2 54 3 98 122 100 123
+228 4 2 54 3 54 99 101 122
+229 4 2 54 3 55 103 99 122
+230 4 2 54 3 49 130 56 105
+231 4 2 54 3 57 132 50 105
+232 4 2 54 3 55 99 98 122
+233 4 2 54 3 103 99 97 53
+234 4 2 54 3 101 15 97 53
+235 4 2 54 3 97 99 101 53
+236 4 2 54 3 105 100 102 123
+237 4 2 54 3 53 97 103 13
+238 4 2 54 3 103 97 99 122
+239 4 2 54 3 53 103 99 13
+240 4 2 54 3 15 101 99 53
+241 4 2 54 3 99 97 101 122
+242 4 2 54 3 105 104 98 123
+243 4 2 54 3 105 57 102 100
+244 4 2 54 3 98 105 56 104
+245 4 2 54 3 12 132 50 57
+246 4 2 54 3 56 6 130 49
+247 4 2 54 3 104 127 23 103
+248 4 2 54 3 129 45 101 102
+249 4 2 54 3 103 55 99 13
+250 4 2 54 3 15 99 101 54
+251 4 2 54 3 14 98 56 104
+252 4 2 54 3 52 101 15 97
+253 4 2 54 3 97 51 103 13
+254 4 2 54 3 105 56 104 123
+255 4 2 54 3 57 105 102 123
+256 4 2 54 3 56 105 98 14
+257 4 2 54 3 100 98 99 122
+258 4 2 54 3 57 12 132 46
+259 4 2 54 3 24 56 6 130
+260 4 2 54 3 105 16 57 100
+261 4 2 54 3 16 102 57 100
+262 4 2 54 3 126 11 48 52
+263 4 2 54 3 5 124 47 51
+264 4 2 54 3 58 98 105 14
+265 4 2 54 3 132 57 46 102
+266 4 2 54 3 24 56 130 104
+267 4 2 54 3 105 98 100 123
+268 4 2 54 3 100 105 16 58
+269 4 2 54 3 11 126 45 101
+270 4 2 54 3 23 124 5 103
+271 4 2 54 3 50 131 49 105
+272 4 2 54 3 45 129 46 102
+273 4 2 54 3 24 127 23 104
+274 4 2 54 3 54 102 16 100
+275 4 2 54 3 104 14 98 55
+276 4 2 54 3 124 5 103 51
+277 4 2 54 3 11 126 101 52
+278 4 2 54 3 105 100 98 58
+279 4 2 54 3 47 125 48 97
+280 4 2 54 3 93 128 129 91
+281 4 2 54 3 93 129 128 123
+282 4 2 54 3 122 129 128 91
+283 4 2 54 3 122 128 129 123
+284 4 2 54 3 122 128 127 90
+285 4 2 54 3 122 127 128 123
+286 4 2 54 3 92 127 128 90
+287 4 2 54 3 92 128 127 123
+288 4 2 54 3 91 125 126 48
+289 4 2 54 3 91 126 125 122
+290 4 2 54 3 97 126 125 48
+291 4 2 54 3 97 125 126 122
+292 4 2 54 3 97 125 124 47
+293 4 2 54 3 97 124 125 122
+294 4 2 54 3 90 124 125 47
+295 4 2 54 3 90 125 124 122
+296 4 2 54 3 105 131 132 50
+297 4 2 54 3 105 132 131 123
+298 4 2 54 3 93 132 131 50
+299 4 2 54 3 93 131 132 123
+300 4 2 54 3 92 131 130 49
+301 4 2 54 3 92 130 131 123
+302 4 2 54 3 105 130 131 49
+303 4 2 54 3 105 131 130 123
+304 5 2 52 1 21 22 59 60 77 79 106 108
+305 5 2 52 1 77 79 106 108 78 80 107 109
+306 5 2 52 1 78 80 107 109 33 34 83 84
+307 5 2 52 1 25 59 22 3 81 106 79 42
+308 5 2 52 1 81 106 79 42 82 107 80 43
+309 5 2 52 1 82 107 80 43 35 83 34 10
+310 5 2 52 1 59 17 18 62 106 67 69 112
+311 5 2 52 1 106 67 69 112 107 68 70 113
+312 5 2 52 1 107 68 70 113 83 28 29 86
+313 5 2 52 1 17 59 25 4 67 106 81 36
+314 5 2 52 1 67 106 81 36 68 107 82 37
+315 5 2 52 1 68 107 82 37 28 83 35 7
+316 5 2 52 1 26 61 62 63 73 110 112 114
+317 5 2 52 1 73 110 112 114 74 111 113 115
+318 5 2 52 1 74 111 113 115 31 85 86 87
+319 5 2 52 1 20 21 60 61 75 77 108 110
+320 5 2 52 1 75 77 108 110 76 78 109 111
+321 5 2 52 1 76 78 109 111 32 33 84 85
+322 5 2 52 1 19 63 62 18 71 114 112 69
+323 5 2 52 1 71 114 112 69 72 115 113 70
+324 5 2 52 1 72 115 113 70 30 87 86 29
+325 5 2 52 1 20 61 26 2 75 110 73 40
+326 5 2 52 1 75 110 73 40 76 111 74 41
+327 5 2 52 1 76 111 74 41 32 85 31 9
+328 5 2 52 1 59 62 61 60 106 112 110 108
+329 5 2 52 1 106 112 110 108 107 113 111 109
+330 5 2 52 1 107 113 111 109 83 86 85 84
+331 5 2 52 1 26 63 19 1 73 114 71 38
+332 5 2 52 1 73 114 71 38 74 115 72 39
+333 5 2 52 1 74 115 72 39 31 87 30 8
+334 6 2 53 2 27 20 2 88 75 40
+335 6 2 53 2 88 75 40 89 76 41
+336 6 2 53 2 89 76 41 44 32 9
+337 6 2 53 2 5 23 27 47 90 88
+338 6 2 53 2 47 90 88 48 91 89
+339 6 2 53 2 48 91 89 11 45 44
+340 6 2 53 2 24 64 23 92 116 90
+341 6 2 53 2 92 116 90 93 117 91
+342 6 2 53 2 93 117 91 46 94 45
+343 6 2 53 2 3 22 65 42 79 118
+344 6 2 53 2 42 79 118 43 80 119
+345 6 2 53 2 43 80 119 10 34 95
+346 6 2 53 2 66 23 64 120 90 116
+347 6 2 53 2 120 90 116 121 91 117
+348 6 2 53 2 121 91 117 96 45 94
+349 6 2 53 2 3 65 6 42 118 49
+350 6 2 53 2 42 118 49 43 119 50
+351 6 2 53 2 43 119 50 10 95 12
+352 6 2 53 2 66 64 21 120 116 77
+353 6 2 53 2 120 116 77 121 117 78
+354 6 2 53 2 121 117 78 96 94 33
+355 6 2 53 2 6 65 24 49 118 92
+356 6 2 53 2 49 118 92 50 119 93
+357 6 2 53 2 50 119 93 12 95 46
+358 6 2 53 2 64 24 65 116 92 118
+359 6 2 53 2 116 92 118 117 93 119
+360 6 2 53 2 117 93 119 94 46 95
+361 6 2 53 2 64 65 22 116 118 79
+362 6 2 53 2 116 118 79 117 119 80
+363 6 2 53 2 117 119 80 94 95 34
+364 6 2 53 2 27 66 20 88 120 75
+365 6 2 53 2 88 120 75 89 121 76
+366 6 2 53 2 89 121 76 44 96 32
+367 6 2 53 2 22 21 64 79 77 116
+368 6 2 53 2 79 77 116 80 78 117
+369 6 2 53 2 80 78 117 34 33 94
+370 6 2 53 2 21 20 66 77 75 120
+371 6 2 53 2 77 75 120 78 76 121
+372 6 2 53 2 78 76 121 33 32 96
+373 6 2 53 2 27 23 66 88 90 120
+374 6 2 53 2 88 90 120 89 91 121
+375 6 2 53 2 89 91 121 44 45 96
+376 7 2 54 3 5 23 90 47 124
+377 7 2 54 3 47 90 91 48 125
+378 7 2 54 3 48 91 45 11 126
+379 7 2 54 3 23 24 92 90 127
+380 7 2 54 3 90 92 93 91 128
+381 7 2 54 3 91 93 46 45 129
+382 7 2 54 3 24 6 49 92 130
+383 7 2 54 3 92 49 50 93 131
+384 7 2 54 3 93 50 12 46 132
+$EndElements
+$Periodic
+9
+0 7 4
+1
+7 4
+0 8 1
+1
+8 1
+1 12 1
+3
+28 17
+29 18
+30 19
+1 13 6
+1
+31 26
+1 15 4
+1
+35 25
+1 35 7
+1
+44 27
+1 37 5
+0
+2 32 1
+5
+86 62
+87 63
+83 59
+84 60
+85 61
+2 54 2
+3
+94 64
+95 65
+96 66
+$EndPeriodic
+)";
+  }
+  return std::dynamic_pointer_cast<const Mesh<Connectivity<3>>>(GmshReader{filename}.mesh());
+}
diff --git a/tests/MeshDataBaseForTests.hpp b/tests/MeshDataBaseForTests.hpp
index 4c34019134f36ce6dbaaffe4528e20c808942926..119348f0e711cb773c46182a46dc835088143a56 100644
--- a/tests/MeshDataBaseForTests.hpp
+++ b/tests/MeshDataBaseForTests.hpp
@@ -23,8 +23,10 @@ class MeshDataBaseForTests
   std::shared_ptr<const Mesh<Connectivity<3>>> m_cartesian_3d_mesh;
 
   std::shared_ptr<const Mesh<Connectivity<2>>> m_hybrid_2d_mesh;
+  std::shared_ptr<const Mesh<Connectivity<3>>> m_hybrid_3d_mesh;
 
   std::shared_ptr<const Mesh<Connectivity<2>>> _buildHybrid2dMesh();
+  std::shared_ptr<const Mesh<Connectivity<3>>> _buildHybrid3dMesh();
 
  public:
   std::shared_ptr<const Mesh<Connectivity<1>>> cartesian1DMesh() const;
@@ -32,6 +34,7 @@ class MeshDataBaseForTests
   std::shared_ptr<const Mesh<Connectivity<3>>> cartesian3DMesh() const;
 
   std::shared_ptr<const Mesh<Connectivity<2>>> hybrid2DMesh() const;
+  std::shared_ptr<const Mesh<Connectivity<3>>> hybrid3DMesh() const;
 
   static const MeshDataBaseForTests& get();
   static void create();
diff --git a/tests/test_DiscreteFunctionInterpoler.cpp b/tests/test_DiscreteFunctionInterpoler.cpp
index ed7ae7c32dcab61467bb0bea2a975a0bd885baaf..c697af1eeaa1c4a39ae3104d226ac173caceb475 100644
--- a/tests/test_DiscreteFunctionInterpoler.cpp
+++ b/tests/test_DiscreteFunctionInterpoler.cpp
@@ -334,302 +334,17 @@ let R3x3_non_linear_1d: R^1 -> R^3x3, x -> (2 * exp(x[0]) * sin(x[0]) + 3, sin(x
   {
     constexpr size_t Dimension = 2;
 
-    SECTION("cartesian grid")
-    {
-      const auto& mesh_2d = MeshDataBaseForTests::get().cartesian2DMesh();
-      auto xj             = MeshDataManager::instance().getMeshData(*mesh_2d).xj();
-
-      std::string_view data = R"(
-import math;
-let B_scalar_non_linear_2d: R^2 -> B, x -> (exp(2 * x[0])< 2*x[1]);
-let N_scalar_non_linear_2d: R^2 -> N, x -> floor(3 * (x[0] + x[1]) * (x[0] + x[1]) + 2);
-let Z_scalar_non_linear_2d: R^2 -> Z, x -> floor(exp(2 * x[0]) - 3 * x[1]);
-let R_scalar_non_linear_2d: R^2 -> R, x -> 2 * exp(x[0]) + 3 * x[1];
-let R1_non_linear_2d: R^2 -> R^1, x -> 2 * exp(x[0]);
-let R2_non_linear_2d: R^2 -> R^2, x -> (2 * exp(x[0]), -3*x[1]);
-let R3_non_linear_2d: R^2 -> R^3, x -> (2 * exp(x[0]) + 3, x[1] - 2, 3);
-let R1x1_non_linear_2d: R^2 -> R^1x1, x -> (2 * exp(x[0]) * sin(x[1]) + 3);
-let R2x2_non_linear_2d: R^2 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3, sin(x[1] - 2 * x[0]), 3, x[1] * x[0]);
-let R3x3_non_linear_2d: R^2 -> R^3x3, x -> (2 * exp(x[0]) * sin(x[1]) + 3, sin(x[1] - 2 * x[0]), 3, x[1] * x[0], -4*x[1], 2*x[0]+1, 3, -6*x[0], exp(x[1]));
-)";
-      TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
-
-      auto ast = ASTBuilder::build(input);
-
-      ASTModulesImporter{*ast};
-      ASTNodeTypeCleaner<language::import_instruction>{*ast};
-
-      ASTSymbolTableBuilder{*ast};
-      ASTNodeDataTypeBuilder{*ast};
-
-      ASTNodeTypeCleaner<language::var_declaration>{*ast};
-      ASTNodeTypeCleaner<language::fct_declaration>{*ast};
-      ASTNodeExpressionBuilder{*ast};
-
-      std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
-
-      TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
-      position.byte = data.size();   // ensure that variables are declared at this point
-
-      SECTION("B_scalar_non_linear_2d")
-      {
-        auto [i_symbol, found] = symbol_table->find("B_scalar_non_linear_2d", position);
-        REQUIRE(found);
-        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
-
-        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
-
-        CellValue<double> cell_value{mesh_2d->connectivity()};
-        parallel_for(
-          cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
-            const TinyVector<Dimension>& x = xj[cell_id];
-            cell_value[cell_id]            = std::exp(2 * x[0]) < 2 * x[1];
-          });
-
-        DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                              function_symbol_id);
-        std::shared_ptr discrete_function = interpoler.interpolate();
-
-        REQUIRE(
-          same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
-      }
-
-      SECTION("N_scalar_non_linear_2d")
-      {
-        auto [i_symbol, found] = symbol_table->find("N_scalar_non_linear_2d", position);
-        REQUIRE(found);
-        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
-
-        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
-
-        CellValue<double> cell_value{mesh_2d->connectivity()};
-        parallel_for(
-          cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
-            const TinyVector<Dimension>& x = xj[cell_id];
-            cell_value[cell_id]            = std::floor(3 * (x[0] + x[1]) * (x[0] + x[1]) + 2);
-          });
-
-        DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                              function_symbol_id);
-        std::shared_ptr discrete_function = interpoler.interpolate();
-
-        REQUIRE(
-          same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
-      }
-
-      SECTION("Z_scalar_non_linear_2d")
-      {
-        auto [i_symbol, found] = symbol_table->find("Z_scalar_non_linear_2d", position);
-        REQUIRE(found);
-        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
-
-        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
-
-        CellValue<double> cell_value{mesh_2d->connectivity()};
-        parallel_for(
-          cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
-            const TinyVector<Dimension>& x = xj[cell_id];
-            cell_value[cell_id]            = std::floor(std::exp(2 * x[0]) - 3 * x[1]);
-          });
-
-        DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                              function_symbol_id);
-        std::shared_ptr discrete_function = interpoler.interpolate();
-
-        REQUIRE(
-          same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
-      }
-
-      SECTION("R_scalar_non_linear_2d")
-      {
-        auto [i_symbol, found] = symbol_table->find("R_scalar_non_linear_2d", position);
-        REQUIRE(found);
-        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
-
-        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
-
-        CellValue<double> cell_value{mesh_2d->connectivity()};
-        parallel_for(
-          cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
-            const TinyVector<Dimension>& x = xj[cell_id];
-            cell_value[cell_id]            = 2 * std::exp(x[0]) + 3 * x[1];
-          });
-
-        DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                              function_symbol_id);
-        std::shared_ptr discrete_function = interpoler.interpolate();
-
-        REQUIRE(
-          same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
-      }
+    std::array mesh_list =   //
+      {std::make_pair(std::string{"cartesian grid"}, MeshDataBaseForTests::get().cartesian2DMesh()),
+       std::make_pair(std::string{"hybrid grid"}, MeshDataBaseForTests::get().hybrid2DMesh())};
 
-      SECTION("R1_non_linear_2d")
+    for (auto [section_name, mesh_3d] : mesh_list) {
+      SECTION(section_name)
       {
-        using DataType = TinyVector<1>;
+        const auto& mesh_2d = MeshDataBaseForTests::get().cartesian2DMesh();
+        auto xj             = MeshDataManager::instance().getMeshData(*mesh_2d).xj();
 
-        auto [i_symbol, found] = symbol_table->find("R1_non_linear_2d", position);
-        REQUIRE(found);
-        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
-
-        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
-
-        CellValue<DataType> cell_value{mesh_2d->connectivity()};
-        parallel_for(
-          cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
-            const TinyVector<Dimension>& x = xj[cell_id];
-            cell_value[cell_id]            = DataType{2 * std::exp(x[0])};
-          });
-
-        DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                              function_symbol_id);
-        std::shared_ptr discrete_function = interpoler.interpolate();
-
-        REQUIRE(same_cell_value(cell_value,
-                                dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
-      }
-
-      SECTION("R2_non_linear_2d")
-      {
-        using DataType = TinyVector<2>;
-
-        auto [i_symbol, found] = symbol_table->find("R2_non_linear_2d", position);
-        REQUIRE(found);
-        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
-
-        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
-
-        CellValue<DataType> cell_value{mesh_2d->connectivity()};
-        parallel_for(
-          cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
-            const TinyVector<Dimension>& x = xj[cell_id];
-            cell_value[cell_id]            = DataType{2 * std::exp(x[0]), -3 * x[1]};
-          });
-
-        DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                              function_symbol_id);
-        std::shared_ptr discrete_function = interpoler.interpolate();
-
-        REQUIRE(same_cell_value(cell_value,
-                                dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
-      }
-
-      SECTION("R3_non_linear_2d")
-      {
-        using DataType = TinyVector<3>;
-
-        auto [i_symbol, found] = symbol_table->find("R3_non_linear_2d", position);
-        REQUIRE(found);
-        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
-
-        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
-
-        CellValue<DataType> cell_value{mesh_2d->connectivity()};
-        parallel_for(
-          cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
-            const TinyVector<Dimension>& x = xj[cell_id];
-            cell_value[cell_id]            = DataType{2 * std::exp(x[0]) + 3, x[1] - 2, 3};
-          });
-
-        DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                              function_symbol_id);
-        std::shared_ptr discrete_function = interpoler.interpolate();
-
-        REQUIRE(same_cell_value(cell_value,
-                                dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
-      }
-
-      SECTION("R1x1_non_linear_2d")
-      {
-        using DataType = TinyMatrix<1>;
-
-        auto [i_symbol, found] = symbol_table->find("R1x1_non_linear_2d", position);
-        REQUIRE(found);
-        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
-
-        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
-
-        CellValue<DataType> cell_value{mesh_2d->connectivity()};
-        parallel_for(
-          cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
-            const TinyVector<Dimension>& x = xj[cell_id];
-            cell_value[cell_id]            = DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3};
-          });
-
-        DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                              function_symbol_id);
-        std::shared_ptr discrete_function = interpoler.interpolate();
-
-        REQUIRE(same_cell_value(cell_value,
-                                dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
-      }
-
-      SECTION("R2x2_non_linear_2d")
-      {
-        using DataType = TinyMatrix<2>;
-
-        auto [i_symbol, found] = symbol_table->find("R2x2_non_linear_2d", position);
-        REQUIRE(found);
-        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
-
-        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
-
-        CellValue<DataType> cell_value{mesh_2d->connectivity()};
-        parallel_for(
-          cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
-            const TinyVector<Dimension>& x = xj[cell_id];
-            cell_value[cell_id] =
-              DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3, std::sin(x[1] - 2 * x[0]), 3, x[1] * x[0]};
-          });
-
-        DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                              function_symbol_id);
-        std::shared_ptr discrete_function = interpoler.interpolate();
-
-        REQUIRE(same_cell_value(cell_value,
-                                dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
-      }
-
-      SECTION("R3x3_non_linear_2d")
-      {
-        using DataType = TinyMatrix<3>;
-
-        auto [i_symbol, found] = symbol_table->find("R3x3_non_linear_2d", position);
-        REQUIRE(found);
-        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
-
-        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
-
-        CellValue<DataType> cell_value{mesh_2d->connectivity()};
-        parallel_for(
-          cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
-            const TinyVector<Dimension>& x = xj[cell_id];
-
-            cell_value[cell_id] = DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3,
-                                           std::sin(x[1] - 2 * x[0]),
-                                           3,
-                                           x[1] * x[0],
-                                           -4 * x[1],
-                                           2 * x[0] + 1,
-                                           3,
-                                           -6 * x[0],
-                                           std::exp(x[1])};
-          });
-
-        DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                              function_symbol_id);
-        std::shared_ptr discrete_function = interpoler.interpolate();
-
-        REQUIRE(same_cell_value(cell_value,
-                                dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
-      }
-    }
-
-    SECTION("hybrid grid")
-    {
-      const auto& mesh_2d = MeshDataBaseForTests::get().hybrid2DMesh();
-      auto xj             = MeshDataManager::instance().getMeshData(*mesh_2d).xj();
-
-      std::string_view data = R"(
+        std::string_view data = R"(
 import math;
 let B_scalar_non_linear_2d: R^2 -> B, x -> (exp(2 * x[0])< 2*x[1]);
 let N_scalar_non_linear_2d: R^2 -> N, x -> floor(3 * (x[0] + x[1]) * (x[0] + x[1]) + 2);
@@ -642,275 +357,276 @@ let R1x1_non_linear_2d: R^2 -> R^1x1, x -> (2 * exp(x[0]) * sin(x[1]) + 3);
 let R2x2_non_linear_2d: R^2 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3, sin(x[1] - 2 * x[0]), 3, x[1] * x[0]);
 let R3x3_non_linear_2d: R^2 -> R^3x3, x -> (2 * exp(x[0]) * sin(x[1]) + 3, sin(x[1] - 2 * x[0]), 3, x[1] * x[0], -4*x[1], 2*x[0]+1, 3, -6*x[0], exp(x[1]));
 )";
-      TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
-
-      auto ast = ASTBuilder::build(input);
-
-      ASTModulesImporter{*ast};
-      ASTNodeTypeCleaner<language::import_instruction>{*ast};
-
-      ASTSymbolTableBuilder{*ast};
-      ASTNodeDataTypeBuilder{*ast};
+        TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+        auto ast = ASTBuilder::build(input);
+
+        ASTModulesImporter{*ast};
+        ASTNodeTypeCleaner<language::import_instruction>{*ast};
 
-      ASTNodeTypeCleaner<language::var_declaration>{*ast};
-      ASTNodeTypeCleaner<language::fct_declaration>{*ast};
-      ASTNodeExpressionBuilder{*ast};
+        ASTSymbolTableBuilder{*ast};
+        ASTNodeDataTypeBuilder{*ast};
+
+        ASTNodeTypeCleaner<language::var_declaration>{*ast};
+        ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+        ASTNodeExpressionBuilder{*ast};
+
+        std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+        TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+        position.byte = data.size();   // ensure that variables are declared at this point
 
-      std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+        SECTION("B_scalar_non_linear_2d")
+        {
+          auto [i_symbol, found] = symbol_table->find("B_scalar_non_linear_2d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          CellValue<double> cell_value{mesh_2d->connectivity()};
+          parallel_for(
+            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+              const TinyVector<Dimension>& x = xj[cell_id];
+              cell_value[cell_id]            = std::exp(2 * x[0]) < 2 * x[1];
+            });
+
+          DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                                function_symbol_id);
+          std::shared_ptr discrete_function = interpoler.interpolate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+        }
+
+        SECTION("N_scalar_non_linear_2d")
+        {
+          auto [i_symbol, found] = symbol_table->find("N_scalar_non_linear_2d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          CellValue<double> cell_value{mesh_2d->connectivity()};
+          parallel_for(
+            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+              const TinyVector<Dimension>& x = xj[cell_id];
+              cell_value[cell_id]            = std::floor(3 * (x[0] + x[1]) * (x[0] + x[1]) + 2);
+            });
+
+          DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                                function_symbol_id);
+          std::shared_ptr discrete_function = interpoler.interpolate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+        }
+
+        SECTION("Z_scalar_non_linear_2d")
+        {
+          auto [i_symbol, found] = symbol_table->find("Z_scalar_non_linear_2d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
 
-      TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
-      position.byte = data.size();   // ensure that variables are declared at this point
-
-      SECTION("B_scalar_non_linear_2d")
-      {
-        auto [i_symbol, found] = symbol_table->find("B_scalar_non_linear_2d", position);
-        REQUIRE(found);
-        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
-
-        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
-
-        CellValue<double> cell_value{mesh_2d->connectivity()};
-        parallel_for(
-          cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
-            const TinyVector<Dimension>& x = xj[cell_id];
-            cell_value[cell_id]            = std::exp(2 * x[0]) < 2 * x[1];
-          });
-
-        DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                              function_symbol_id);
-        std::shared_ptr discrete_function = interpoler.interpolate();
-
-        REQUIRE(
-          same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
-      }
-
-      SECTION("N_scalar_non_linear_2d")
-      {
-        auto [i_symbol, found] = symbol_table->find("N_scalar_non_linear_2d", position);
-        REQUIRE(found);
-        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
-
-        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
-
-        CellValue<double> cell_value{mesh_2d->connectivity()};
-        parallel_for(
-          cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
-            const TinyVector<Dimension>& x = xj[cell_id];
-            cell_value[cell_id]            = std::floor(3 * (x[0] + x[1]) * (x[0] + x[1]) + 2);
-          });
-
-        DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                              function_symbol_id);
-        std::shared_ptr discrete_function = interpoler.interpolate();
-
-        REQUIRE(
-          same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
-      }
-
-      SECTION("Z_scalar_non_linear_2d")
-      {
-        auto [i_symbol, found] = symbol_table->find("Z_scalar_non_linear_2d", position);
-        REQUIRE(found);
-        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
-
-        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
-
-        CellValue<double> cell_value{mesh_2d->connectivity()};
-        parallel_for(
-          cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
-            const TinyVector<Dimension>& x = xj[cell_id];
-            cell_value[cell_id]            = std::floor(std::exp(2 * x[0]) - 3 * x[1]);
-          });
-
-        DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                              function_symbol_id);
-        std::shared_ptr discrete_function = interpoler.interpolate();
-
-        REQUIRE(
-          same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
-      }
-
-      SECTION("R_scalar_non_linear_2d")
-      {
-        auto [i_symbol, found] = symbol_table->find("R_scalar_non_linear_2d", position);
-        REQUIRE(found);
-        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
-
-        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
-
-        CellValue<double> cell_value{mesh_2d->connectivity()};
-        parallel_for(
-          cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
-            const TinyVector<Dimension>& x = xj[cell_id];
-            cell_value[cell_id]            = 2 * std::exp(x[0]) + 3 * x[1];
-          });
-
-        DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                              function_symbol_id);
-        std::shared_ptr discrete_function = interpoler.interpolate();
-
-        REQUIRE(
-          same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
-      }
-
-      SECTION("R1_non_linear_2d")
-      {
-        using DataType = TinyVector<1>;
-
-        auto [i_symbol, found] = symbol_table->find("R1_non_linear_2d", position);
-        REQUIRE(found);
-        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
-
-        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
-
-        CellValue<DataType> cell_value{mesh_2d->connectivity()};
-        parallel_for(
-          cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
-            const TinyVector<Dimension>& x = xj[cell_id];
-            cell_value[cell_id]            = DataType{2 * std::exp(x[0])};
-          });
-
-        DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                              function_symbol_id);
-        std::shared_ptr discrete_function = interpoler.interpolate();
-
-        REQUIRE(same_cell_value(cell_value,
-                                dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
-      }
-
-      SECTION("R2_non_linear_2d")
-      {
-        using DataType = TinyVector<2>;
-
-        auto [i_symbol, found] = symbol_table->find("R2_non_linear_2d", position);
-        REQUIRE(found);
-        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
-
-        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
-
-        CellValue<DataType> cell_value{mesh_2d->connectivity()};
-        parallel_for(
-          cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
-            const TinyVector<Dimension>& x = xj[cell_id];
-            cell_value[cell_id]            = DataType{2 * std::exp(x[0]), -3 * x[1]};
-          });
-
-        DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                              function_symbol_id);
-        std::shared_ptr discrete_function = interpoler.interpolate();
-
-        REQUIRE(same_cell_value(cell_value,
-                                dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
-      }
-
-      SECTION("R3_non_linear_2d")
-      {
-        using DataType = TinyVector<3>;
-
-        auto [i_symbol, found] = symbol_table->find("R3_non_linear_2d", position);
-        REQUIRE(found);
-        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
-
-        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
-
-        CellValue<DataType> cell_value{mesh_2d->connectivity()};
-        parallel_for(
-          cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
-            const TinyVector<Dimension>& x = xj[cell_id];
-            cell_value[cell_id]            = DataType{2 * std::exp(x[0]) + 3, x[1] - 2, 3};
-          });
-
-        DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                              function_symbol_id);
-        std::shared_ptr discrete_function = interpoler.interpolate();
-
-        REQUIRE(same_cell_value(cell_value,
-                                dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
-      }
-
-      SECTION("R1x1_non_linear_2d")
-      {
-        using DataType = TinyMatrix<1>;
-
-        auto [i_symbol, found] = symbol_table->find("R1x1_non_linear_2d", position);
-        REQUIRE(found);
-        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
-
-        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
-
-        CellValue<DataType> cell_value{mesh_2d->connectivity()};
-        parallel_for(
-          cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
-            const TinyVector<Dimension>& x = xj[cell_id];
-            cell_value[cell_id]            = DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3};
-          });
-
-        DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                              function_symbol_id);
-        std::shared_ptr discrete_function = interpoler.interpolate();
-
-        REQUIRE(same_cell_value(cell_value,
-                                dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
-      }
-
-      SECTION("R2x2_non_linear_2d")
-      {
-        using DataType = TinyMatrix<2>;
-
-        auto [i_symbol, found] = symbol_table->find("R2x2_non_linear_2d", position);
-        REQUIRE(found);
-        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
-
-        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
-
-        CellValue<DataType> cell_value{mesh_2d->connectivity()};
-        parallel_for(
-          cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
-            const TinyVector<Dimension>& x = xj[cell_id];
-            cell_value[cell_id] =
-              DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3, std::sin(x[1] - 2 * x[0]), 3, x[1] * x[0]};
-          });
-
-        DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                              function_symbol_id);
-        std::shared_ptr discrete_function = interpoler.interpolate();
-
-        REQUIRE(same_cell_value(cell_value,
-                                dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
-      }
-
-      SECTION("R3x3_non_linear_2d")
-      {
-        using DataType = TinyMatrix<3>;
-
-        auto [i_symbol, found] = symbol_table->find("R3x3_non_linear_2d", position);
-        REQUIRE(found);
-        REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
-
-        FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
-
-        CellValue<DataType> cell_value{mesh_2d->connectivity()};
-        parallel_for(
-          cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
-            const TinyVector<Dimension>& x = xj[cell_id];
-
-            cell_value[cell_id] = DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3,
-                                           std::sin(x[1] - 2 * x[0]),
-                                           3,
-                                           x[1] * x[0],
-                                           -4 * x[1],
-                                           2 * x[0] + 1,
-                                           3,
-                                           -6 * x[0],
-                                           std::exp(x[1])};
-          });
-
-        DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                              function_symbol_id);
-        std::shared_ptr discrete_function = interpoler.interpolate();
-
-        REQUIRE(same_cell_value(cell_value,
-                                dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          CellValue<double> cell_value{mesh_2d->connectivity()};
+          parallel_for(
+            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+              const TinyVector<Dimension>& x = xj[cell_id];
+              cell_value[cell_id]            = std::floor(std::exp(2 * x[0]) - 3 * x[1]);
+            });
+
+          DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                                function_symbol_id);
+          std::shared_ptr discrete_function = interpoler.interpolate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+        }
+
+        SECTION("R_scalar_non_linear_2d")
+        {
+          auto [i_symbol, found] = symbol_table->find("R_scalar_non_linear_2d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          CellValue<double> cell_value{mesh_2d->connectivity()};
+          parallel_for(
+            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+              const TinyVector<Dimension>& x = xj[cell_id];
+              cell_value[cell_id]            = 2 * std::exp(x[0]) + 3 * x[1];
+            });
+
+          DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                                function_symbol_id);
+          std::shared_ptr discrete_function = interpoler.interpolate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+        }
+
+        SECTION("R1_non_linear_2d")
+        {
+          using DataType = TinyVector<1>;
+
+          auto [i_symbol, found] = symbol_table->find("R1_non_linear_2d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          CellValue<DataType> cell_value{mesh_2d->connectivity()};
+          parallel_for(
+            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+              const TinyVector<Dimension>& x = xj[cell_id];
+              cell_value[cell_id]            = DataType{2 * std::exp(x[0])};
+            });
+
+          DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                                function_symbol_id);
+          std::shared_ptr discrete_function = interpoler.interpolate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+        }
+
+        SECTION("R2_non_linear_2d")
+        {
+          using DataType = TinyVector<2>;
+
+          auto [i_symbol, found] = symbol_table->find("R2_non_linear_2d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          CellValue<DataType> cell_value{mesh_2d->connectivity()};
+          parallel_for(
+            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+              const TinyVector<Dimension>& x = xj[cell_id];
+              cell_value[cell_id]            = DataType{2 * std::exp(x[0]), -3 * x[1]};
+            });
+
+          DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                                function_symbol_id);
+          std::shared_ptr discrete_function = interpoler.interpolate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+        }
+
+        SECTION("R3_non_linear_2d")
+        {
+          using DataType = TinyVector<3>;
+
+          auto [i_symbol, found] = symbol_table->find("R3_non_linear_2d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          CellValue<DataType> cell_value{mesh_2d->connectivity()};
+          parallel_for(
+            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+              const TinyVector<Dimension>& x = xj[cell_id];
+              cell_value[cell_id]            = DataType{2 * std::exp(x[0]) + 3, x[1] - 2, 3};
+            });
+
+          DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                                function_symbol_id);
+          std::shared_ptr discrete_function = interpoler.interpolate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+        }
+
+        SECTION("R1x1_non_linear_2d")
+        {
+          using DataType = TinyMatrix<1>;
+
+          auto [i_symbol, found] = symbol_table->find("R1x1_non_linear_2d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          CellValue<DataType> cell_value{mesh_2d->connectivity()};
+          parallel_for(
+            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+              const TinyVector<Dimension>& x = xj[cell_id];
+              cell_value[cell_id]            = DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3};
+            });
+
+          DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                                function_symbol_id);
+          std::shared_ptr discrete_function = interpoler.interpolate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+        }
+
+        SECTION("R2x2_non_linear_2d")
+        {
+          using DataType = TinyMatrix<2>;
+
+          auto [i_symbol, found] = symbol_table->find("R2x2_non_linear_2d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          CellValue<DataType> cell_value{mesh_2d->connectivity()};
+          parallel_for(
+            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+              const TinyVector<Dimension>& x = xj[cell_id];
+              cell_value[cell_id] =
+                DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3, std::sin(x[1] - 2 * x[0]), 3, x[1] * x[0]};
+            });
+
+          DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                                function_symbol_id);
+          std::shared_ptr discrete_function = interpoler.interpolate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+        }
+
+        SECTION("R3x3_non_linear_2d")
+        {
+          using DataType = TinyMatrix<3>;
+
+          auto [i_symbol, found] = symbol_table->find("R3x3_non_linear_2d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          CellValue<DataType> cell_value{mesh_2d->connectivity()};
+          parallel_for(
+            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+              const TinyVector<Dimension>& x = xj[cell_id];
+
+              cell_value[cell_id] = DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3,
+                                             std::sin(x[1] - 2 * x[0]),
+                                             3,
+                                             x[1] * x[0],
+                                             -4 * x[1],
+                                             2 * x[0] + 1,
+                                             3,
+                                             -6 * x[0],
+                                             std::exp(x[1])};
+            });
+
+          DiscreteFunctionInterpoler interpoler(mesh_2d, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                                function_symbol_id);
+          std::shared_ptr discrete_function = interpoler.interpolate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+        }
       }
     }
   }
@@ -919,10 +635,16 @@ let R3x3_non_linear_2d: R^2 -> R^3x3, x -> (2 * exp(x[0]) * sin(x[1]) + 3, sin(x
   {
     constexpr size_t Dimension = 3;
 
-    const auto& mesh_3d = MeshDataBaseForTests::get().cartesian3DMesh();
-    auto xj             = MeshDataManager::instance().getMeshData(*mesh_3d).xj();
+    std::array mesh_list =   //
+      {std::make_pair(std::string{"cartesian grid"}, MeshDataBaseForTests::get().cartesian3DMesh()),
+       std::make_pair(std::string{"hybrid grid"}, MeshDataBaseForTests::get().hybrid3DMesh())};
 
-    std::string_view data = R"(
+    for (auto [section_name, mesh_3d] : mesh_list) {
+      SECTION(section_name)
+      {
+        auto xj = MeshDataManager::instance().getMeshData(*mesh_3d).xj();
+
+        std::string_view data = R"(
 import math;
 let B_scalar_non_linear_3d: R^3 -> B, x -> (exp(2 * x[0])< 2*x[1]+x[2]);
 let N_scalar_non_linear_3d: R^3 -> N, x -> floor(3 * (x[0] + x[1]) * (x[0] + x[1]) + x[2] * x[2]);
@@ -935,275 +657,277 @@ let R1x1_non_linear_3d: R^3 -> R^1x1, x -> (2 * exp(x[0]) * sin(x[1]) + 3 * x[2]
 let R2x2_non_linear_3d: R^3 -> R^2x2, x -> (2 * exp(x[0]) * sin(x[1]) + 3, sin(x[2] - 2 * x[0]), 3, x[1] * x[0] - x[2]);
 let R3x3_non_linear_3d: R^3 -> R^3x3, x -> (2 * exp(x[0]) * sin(x[1]) + 3, sin(x[1] - 2 * x[2]), 3, x[1] * x[2], -4*x[1], 2*x[2]+1, 3, -6*x[2], exp(x[1] + x[2]));
 )";
-    TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
-
-    auto ast = ASTBuilder::build(input);
-
-    ASTModulesImporter{*ast};
-    ASTNodeTypeCleaner<language::import_instruction>{*ast};
+        TAO_PEGTL_NAMESPACE::string_input input{data, "test.pgs"};
+
+        auto ast = ASTBuilder::build(input);
+
+        ASTModulesImporter{*ast};
+        ASTNodeTypeCleaner<language::import_instruction>{*ast};
 
-    ASTSymbolTableBuilder{*ast};
-    ASTNodeDataTypeBuilder{*ast};
+        ASTSymbolTableBuilder{*ast};
+        ASTNodeDataTypeBuilder{*ast};
+
+        ASTNodeTypeCleaner<language::var_declaration>{*ast};
+        ASTNodeTypeCleaner<language::fct_declaration>{*ast};
+        ASTNodeExpressionBuilder{*ast};
+
+        std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
+
+        TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
+        position.byte = data.size();   // ensure that variables are declared at this point
 
-    ASTNodeTypeCleaner<language::var_declaration>{*ast};
-    ASTNodeTypeCleaner<language::fct_declaration>{*ast};
-    ASTNodeExpressionBuilder{*ast};
-
-    std::shared_ptr<SymbolTable> symbol_table = ast->m_symbol_table;
-
-    TAO_PEGTL_NAMESPACE::position position{TAO_PEGTL_NAMESPACE::internal::iterator{"fixture"}, "fixture"};
-    position.byte = data.size();   // ensure that variables are declared at this point
-
-    SECTION("B_scalar_non_linear_3d")
-    {
-      auto [i_symbol, found] = symbol_table->find("B_scalar_non_linear_3d", position);
-      REQUIRE(found);
-      REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
-
-      FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+        SECTION("B_scalar_non_linear_3d")
+        {
+          auto [i_symbol, found] = symbol_table->find("B_scalar_non_linear_3d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          CellValue<double> cell_value{mesh_3d->connectivity()};
+          parallel_for(
+            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+              const TinyVector<Dimension>& x = xj[cell_id];
+              cell_value[cell_id]            = std::exp(2 * x[0]) < 2 * x[1] + x[2];
+            });
+
+          DiscreteFunctionInterpoler interpoler(mesh_3d, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                                function_symbol_id);
+          std::shared_ptr discrete_function = interpoler.interpolate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+        }
+
+        SECTION("N_scalar_non_linear_3d")
+        {
+          auto [i_symbol, found] = symbol_table->find("N_scalar_non_linear_3d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          CellValue<double> cell_value{mesh_3d->connectivity()};
+          parallel_for(
+            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+              const TinyVector<Dimension>& x = xj[cell_id];
+              cell_value[cell_id]            = std::floor(3 * (x[0] + x[1]) * (x[0] + x[1]) + x[2] * x[2]);
+            });
+
+          DiscreteFunctionInterpoler interpoler(mesh_3d, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                                function_symbol_id);
+          std::shared_ptr discrete_function = interpoler.interpolate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+        }
+
+        SECTION("Z_scalar_non_linear_3d")
+        {
+          auto [i_symbol, found] = symbol_table->find("Z_scalar_non_linear_3d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
 
-      CellValue<double> cell_value{mesh_3d->connectivity()};
-      parallel_for(
-        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
-          const TinyVector<Dimension>& x = xj[cell_id];
-          cell_value[cell_id]            = std::exp(2 * x[0]) < 2 * x[1] + x[2];
-        });
-
-      DiscreteFunctionInterpoler interpoler(mesh_3d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                            function_symbol_id);
-      std::shared_ptr discrete_function = interpoler.interpolate();
-
-      REQUIRE(
-        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
-    }
-
-    SECTION("N_scalar_non_linear_3d")
-    {
-      auto [i_symbol, found] = symbol_table->find("N_scalar_non_linear_3d", position);
-      REQUIRE(found);
-      REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
-
-      FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
-
-      CellValue<double> cell_value{mesh_3d->connectivity()};
-      parallel_for(
-        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
-          const TinyVector<Dimension>& x = xj[cell_id];
-          cell_value[cell_id]            = std::floor(3 * (x[0] + x[1]) * (x[0] + x[1]) + x[2] * x[2]);
-        });
-
-      DiscreteFunctionInterpoler interpoler(mesh_3d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                            function_symbol_id);
-      std::shared_ptr discrete_function = interpoler.interpolate();
-
-      REQUIRE(
-        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
-    }
-
-    SECTION("Z_scalar_non_linear_3d")
-    {
-      auto [i_symbol, found] = symbol_table->find("Z_scalar_non_linear_3d", position);
-      REQUIRE(found);
-      REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
-
-      FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
-
-      CellValue<double> cell_value{mesh_3d->connectivity()};
-      parallel_for(
-        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
-          const TinyVector<Dimension>& x = xj[cell_id];
-          cell_value[cell_id]            = std::floor(std::exp(2 * x[0]) - 3 * x[1] + x[2]);
-        });
-
-      DiscreteFunctionInterpoler interpoler(mesh_3d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                            function_symbol_id);
-      std::shared_ptr discrete_function = interpoler.interpolate();
-
-      REQUIRE(
-        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
-    }
-
-    SECTION("R_scalar_non_linear_3d")
-    {
-      auto [i_symbol, found] = symbol_table->find("R_scalar_non_linear_3d", position);
-      REQUIRE(found);
-      REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
-
-      FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
-
-      CellValue<double> cell_value{mesh_3d->connectivity()};
-      parallel_for(
-        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
-          const TinyVector<Dimension>& x = xj[cell_id];
-          cell_value[cell_id]            = 2 * std::exp(x[0] + x[2]) + 3 * x[1];
-        });
-
-      DiscreteFunctionInterpoler interpoler(mesh_3d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                            function_symbol_id);
-      std::shared_ptr discrete_function = interpoler.interpolate();
-
-      REQUIRE(
-        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
-    }
-
-    SECTION("R1_non_linear_3d")
-    {
-      using DataType = TinyVector<1>;
-
-      auto [i_symbol, found] = symbol_table->find("R1_non_linear_3d", position);
-      REQUIRE(found);
-      REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
-
-      FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
-
-      CellValue<DataType> cell_value{mesh_3d->connectivity()};
-      parallel_for(
-        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
-          const TinyVector<Dimension>& x = xj[cell_id];
-          cell_value[cell_id]            = DataType{2 * std::exp(x[0]) + std::sin(x[1] + x[2])};
-        });
-
-      DiscreteFunctionInterpoler interpoler(mesh_3d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                            function_symbol_id);
-      std::shared_ptr discrete_function = interpoler.interpolate();
-
-      REQUIRE(
-        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
-    }
-
-    SECTION("R2_non_linear_3d")
-    {
-      using DataType = TinyVector<2>;
-
-      auto [i_symbol, found] = symbol_table->find("R2_non_linear_3d", position);
-      REQUIRE(found);
-      REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
-
-      FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
-
-      CellValue<DataType> cell_value{mesh_3d->connectivity()};
-      parallel_for(
-        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
-          const TinyVector<Dimension>& x = xj[cell_id];
-          cell_value[cell_id]            = DataType{2 * std::exp(x[0]), -3 * x[1] * x[2]};
-        });
-
-      DiscreteFunctionInterpoler interpoler(mesh_3d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                            function_symbol_id);
-      std::shared_ptr discrete_function = interpoler.interpolate();
-
-      REQUIRE(
-        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
-    }
-
-    SECTION("R3_non_linear_3d")
-    {
-      using DataType = TinyVector<3>;
-
-      auto [i_symbol, found] = symbol_table->find("R3_non_linear_3d", position);
-      REQUIRE(found);
-      REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
-
-      FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
-
-      CellValue<DataType> cell_value{mesh_3d->connectivity()};
-      parallel_for(
-        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
-          const TinyVector<Dimension>& x = xj[cell_id];
-          cell_value[cell_id]            = DataType{2 * std::exp(x[0]) + 3, x[1] - 2, 3 * x[2]};
-        });
-
-      DiscreteFunctionInterpoler interpoler(mesh_3d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                            function_symbol_id);
-      std::shared_ptr discrete_function = interpoler.interpolate();
-
-      REQUIRE(
-        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
-    }
-
-    SECTION("R1x1_non_linear_3d")
-    {
-      using DataType = TinyMatrix<1>;
-
-      auto [i_symbol, found] = symbol_table->find("R1x1_non_linear_3d", position);
-      REQUIRE(found);
-      REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
-
-      FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
-
-      CellValue<DataType> cell_value{mesh_3d->connectivity()};
-      parallel_for(
-        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
-          const TinyVector<Dimension>& x = xj[cell_id];
-          cell_value[cell_id]            = DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3 * x[2]};
-        });
-
-      DiscreteFunctionInterpoler interpoler(mesh_3d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                            function_symbol_id);
-      std::shared_ptr discrete_function = interpoler.interpolate();
-
-      REQUIRE(
-        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
-    }
-
-    SECTION("R2x2_non_linear_3d")
-    {
-      using DataType = TinyMatrix<2>;
-
-      auto [i_symbol, found] = symbol_table->find("R2x2_non_linear_3d", position);
-      REQUIRE(found);
-      REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
-
-      FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
-
-      CellValue<DataType> cell_value{mesh_3d->connectivity()};
-      parallel_for(
-        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
-          const TinyVector<Dimension>& x = xj[cell_id];
-          cell_value[cell_id] =
-            DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3, std::sin(x[2] - 2 * x[0]), 3, x[1] * x[0] - x[2]};
-        });
-
-      DiscreteFunctionInterpoler interpoler(mesh_3d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                            function_symbol_id);
-      std::shared_ptr discrete_function = interpoler.interpolate();
-
-      REQUIRE(
-        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
-    }
-
-    SECTION("R3x3_non_linear_3d")
-    {
-      using DataType = TinyMatrix<3>;
-
-      auto [i_symbol, found] = symbol_table->find("R3x3_non_linear_3d", position);
-      REQUIRE(found);
-      REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
-
-      FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
-
-      CellValue<DataType> cell_value{mesh_3d->connectivity()};
-      parallel_for(
-        cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
-          const TinyVector<Dimension>& x = xj[cell_id];
-
-          cell_value[cell_id] = DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3,
-                                         std::sin(x[1] - 2 * x[2]),
-                                         3,
-                                         x[1] * x[2],
-                                         -4 * x[1],
-                                         2 * x[2] + 1,
-                                         3,
-                                         -6 * x[2],
-                                         std::exp(x[1] + x[2])};
-        });
-
-      DiscreteFunctionInterpoler interpoler(mesh_3d, std::make_shared<DiscreteFunctionDescriptorP0>(),
-                                            function_symbol_id);
-      std::shared_ptr discrete_function = interpoler.interpolate();
-
-      REQUIRE(
-        same_cell_value(cell_value, dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          CellValue<double> cell_value{mesh_3d->connectivity()};
+          parallel_for(
+            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+              const TinyVector<Dimension>& x = xj[cell_id];
+              cell_value[cell_id]            = std::floor(std::exp(2 * x[0]) - 3 * x[1] + x[2]);
+            });
+
+          DiscreteFunctionInterpoler interpoler(mesh_3d, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                                function_symbol_id);
+          std::shared_ptr discrete_function = interpoler.interpolate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+        }
+
+        SECTION("R_scalar_non_linear_3d")
+        {
+          auto [i_symbol, found] = symbol_table->find("R_scalar_non_linear_3d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          CellValue<double> cell_value{mesh_3d->connectivity()};
+          parallel_for(
+            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+              const TinyVector<Dimension>& x = xj[cell_id];
+              cell_value[cell_id]            = 2 * std::exp(x[0] + x[2]) + 3 * x[1];
+            });
+
+          DiscreteFunctionInterpoler interpoler(mesh_3d, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                                function_symbol_id);
+          std::shared_ptr discrete_function = interpoler.interpolate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*discrete_function)));
+        }
+
+        SECTION("R1_non_linear_3d")
+        {
+          using DataType = TinyVector<1>;
+
+          auto [i_symbol, found] = symbol_table->find("R1_non_linear_3d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          CellValue<DataType> cell_value{mesh_3d->connectivity()};
+          parallel_for(
+            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+              const TinyVector<Dimension>& x = xj[cell_id];
+              cell_value[cell_id]            = DataType{2 * std::exp(x[0]) + std::sin(x[1] + x[2])};
+            });
+
+          DiscreteFunctionInterpoler interpoler(mesh_3d, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                                function_symbol_id);
+          std::shared_ptr discrete_function = interpoler.interpolate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+        }
+
+        SECTION("R2_non_linear_3d")
+        {
+          using DataType = TinyVector<2>;
+
+          auto [i_symbol, found] = symbol_table->find("R2_non_linear_3d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          CellValue<DataType> cell_value{mesh_3d->connectivity()};
+          parallel_for(
+            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+              const TinyVector<Dimension>& x = xj[cell_id];
+              cell_value[cell_id]            = DataType{2 * std::exp(x[0]), -3 * x[1] * x[2]};
+            });
+
+          DiscreteFunctionInterpoler interpoler(mesh_3d, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                                function_symbol_id);
+          std::shared_ptr discrete_function = interpoler.interpolate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+        }
+
+        SECTION("R3_non_linear_3d")
+        {
+          using DataType = TinyVector<3>;
+
+          auto [i_symbol, found] = symbol_table->find("R3_non_linear_3d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          CellValue<DataType> cell_value{mesh_3d->connectivity()};
+          parallel_for(
+            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+              const TinyVector<Dimension>& x = xj[cell_id];
+              cell_value[cell_id]            = DataType{2 * std::exp(x[0]) + 3, x[1] - 2, 3 * x[2]};
+            });
+
+          DiscreteFunctionInterpoler interpoler(mesh_3d, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                                function_symbol_id);
+          std::shared_ptr discrete_function = interpoler.interpolate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+        }
+
+        SECTION("R1x1_non_linear_3d")
+        {
+          using DataType = TinyMatrix<1>;
+
+          auto [i_symbol, found] = symbol_table->find("R1x1_non_linear_3d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          CellValue<DataType> cell_value{mesh_3d->connectivity()};
+          parallel_for(
+            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+              const TinyVector<Dimension>& x = xj[cell_id];
+              cell_value[cell_id]            = DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3 * x[2]};
+            });
+
+          DiscreteFunctionInterpoler interpoler(mesh_3d, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                                function_symbol_id);
+          std::shared_ptr discrete_function = interpoler.interpolate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+        }
+
+        SECTION("R2x2_non_linear_3d")
+        {
+          using DataType = TinyMatrix<2>;
+
+          auto [i_symbol, found] = symbol_table->find("R2x2_non_linear_3d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          CellValue<DataType> cell_value{mesh_3d->connectivity()};
+          parallel_for(
+            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+              const TinyVector<Dimension>& x = xj[cell_id];
+              cell_value[cell_id] =
+                DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3, std::sin(x[2] - 2 * x[0]), 3, x[1] * x[0] - x[2]};
+            });
+
+          DiscreteFunctionInterpoler interpoler(mesh_3d, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                                function_symbol_id);
+          std::shared_ptr discrete_function = interpoler.interpolate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+        }
+
+        SECTION("R3x3_non_linear_3d")
+        {
+          using DataType = TinyMatrix<3>;
+
+          auto [i_symbol, found] = symbol_table->find("R3x3_non_linear_3d", position);
+          REQUIRE(found);
+          REQUIRE(i_symbol->attributes().dataType() == ASTNodeDataType::function_t);
+
+          FunctionSymbolId function_symbol_id(std::get<uint64_t>(i_symbol->attributes().value()), symbol_table);
+
+          CellValue<DataType> cell_value{mesh_3d->connectivity()};
+          parallel_for(
+            cell_value.numberOfItems(), PUGS_LAMBDA(const CellId cell_id) {
+              const TinyVector<Dimension>& x = xj[cell_id];
+
+              cell_value[cell_id] = DataType{2 * std::exp(x[0]) * std::sin(x[1]) + 3,
+                                             std::sin(x[1] - 2 * x[2]),
+                                             3,
+                                             x[1] * x[2],
+                                             -4 * x[1],
+                                             2 * x[2] + 1,
+                                             3,
+                                             -6 * x[2],
+                                             std::exp(x[1] + x[2])};
+            });
+
+          DiscreteFunctionInterpoler interpoler(mesh_3d, std::make_shared<DiscreteFunctionDescriptorP0>(),
+                                                function_symbol_id);
+          std::shared_ptr discrete_function = interpoler.interpolate();
+
+          REQUIRE(same_cell_value(cell_value,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, DataType>&>(*discrete_function)));
+        }
+      }
     }
   }
 }