1 // Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced 2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files 3 // LICENSE and NOTICE for details. LLNL-CODE-806117. 4 // 5 // This file is part of the MFEM library. For more information and source code 6 // availability visit https://mfem.org. 7 // 8 // MFEM is free software; you can redistribute it and/or modify it under the 9 // terms of the BSD-3 license. We welcome feedback and contributions, see file 10 // CONTRIBUTING.md for details. 11 12 #ifndef MFEM_TEMPLATE_LAYOUT 13 #define MFEM_TEMPLATE_LAYOUT 14 15 #include "../config/tconfig.hpp" 16 #include "../fem/fespace.hpp" 17 #include "../general/backends.hpp" 18 19 namespace mfem 20 { 21 22 // Layout classes 23 24 template <int N1, int S1> 25 struct OffsetStridedLayout1D; 26 template <int N1, int S1, int N2, int S2> 27 struct StridedLayout2D; 28 29 template <int N1, int S1> 30 struct StridedLayout1D 31 { 32 static const int rank = 1; 33 static const int dim_1 = N1; 34 static const int size = N1; 35 indmfem::StridedLayout1D36 MFEM_HOST_DEVICE static inline int ind(int i1) 37 { 38 return S1*i1; 39 } 40 41 template <int M1> submfem::StridedLayout1D42 static OffsetStridedLayout1D<M1,S1> sub(int o1) 43 { 44 return OffsetStridedLayout1D<M1,S1>(S1*o1); 45 } 46 47 // reshape methods 48 49 template <int N1_1, int N1_2> split_1mfem::StridedLayout1D50 static StridedLayout2D<N1_1,S1,N1_2,S1*N1_1> split_1() 51 { 52 // S1*i1 == S1*(i1_1+N1_1*i1_2) 53 MFEM_STATIC_ASSERT(N1_1*N1_2 == N1, "invalid dimensions"); 54 return StridedLayout2D<N1_1,S1,N1_2,S1*N1_1>(); 55 } 56 }; 57 58 template <int N1, int S1, int N2, int S2> 59 struct OffsetStridedLayout2D; 60 61 template <int N1, int S1> 62 struct OffsetStridedLayout1D 63 { 64 static const int rank = 1; 65 static const int dim_1 = N1; 66 static const int size = N1; 67 68 int offset; 69 OffsetStridedLayout1Dmfem::OffsetStridedLayout1D70 OffsetStridedLayout1D() { } OffsetStridedLayout1Dmfem::OffsetStridedLayout1D71 OffsetStridedLayout1D(int offset_) : offset(offset_) { } indmfem::OffsetStridedLayout1D72 MFEM_HOST_DEVICE inline int ind(int i1) const 73 { 74 return offset+S1*i1; 75 } 76 77 template <int M1> submfem::OffsetStridedLayout1D78 OffsetStridedLayout1D<M1,S1> sub(int o1) const 79 { 80 return OffsetStridedLayout1D<M1,S1>(offset+S1*o1); 81 } 82 83 // reshape methods 84 85 template <int N1_1, int N1_2> split_1mfem::OffsetStridedLayout1D86 OffsetStridedLayout2D<N1_1,S1,N1_2,S1*N1_1> split_1() const 87 { 88 // S1*i1 == S1*(i1_1+N1_1*i1_2) 89 MFEM_STATIC_ASSERT(N1_1*N1_2 == N1, "invalid dimensions"); 90 return OffsetStridedLayout2D<N1_1,S1,N1_2,S1*N1_1>(offset); 91 } 92 }; 93 94 template <int N1, int S1, int N2, int S2, int N3, int S3> 95 struct StridedLayout3D; 96 template <int N1, int S1, int N2, int S2, int N3, int S3, int N4, int S4> 97 struct StridedLayout4D; 98 99 template <int N1, int S1, int N2, int S2> 100 struct StridedLayout2D 101 { 102 static const int rank = 2; 103 static const int dim_1 = N1; 104 static const int dim_2 = N2; 105 static const int size = N1*N2; 106 indmfem::StridedLayout2D107 MFEM_HOST_DEVICE static inline int ind(int i1, int i2) 108 { 109 return (S1*i1+S2*i2); 110 } ind1mfem::StridedLayout2D111 static OffsetStridedLayout1D<N2,S2> ind1(int i1) 112 { 113 return OffsetStridedLayout1D<N2,S2>(S1*i1); 114 } ind2mfem::StridedLayout2D115 static OffsetStridedLayout1D<N1,S1> ind2(int i2) 116 { 117 return OffsetStridedLayout1D<N1,S1>(S2*i2); 118 } 119 120 template <int M1, int M2> submfem::StridedLayout2D121 static OffsetStridedLayout2D<M1,S1,M2,S2> sub(int o1, int o2) 122 { 123 return OffsetStridedLayout2D<M1,S1,M2,S2>(S1*o1+S2*o2); 124 } 125 126 // reshape methods 127 128 template <int N1_1, int N1_2> split_1mfem::StridedLayout2D129 static StridedLayout3D<N1_1,S1,N1_2,S1*N1_1,N2,S2> split_1() 130 { 131 // S1*i1+S2*i2 == S1*(i1_1+N1_1*i1_2)+S2*i2 132 MFEM_STATIC_ASSERT(N1_1*N1_2 == N1, "invalid dimensions"); 133 return StridedLayout3D<N1_1,S1,N1_2,S1*N1_1,N2,S2>(); 134 } 135 template <int N2_1, int N2_2> split_2mfem::StridedLayout2D136 static StridedLayout3D<N1,S1,N2_1,S2,N2_2,S2*N2_1> split_2() 137 { 138 // S1*i1+S2*i2 == S1*i1+S2*(i2_1*N2_1*i2_2) 139 MFEM_STATIC_ASSERT(N2_1*N2_2 == N2, "invalid dimensions"); 140 return StridedLayout3D<N1,S1,N2_1,S2,N2_2,S2*N2_1>(); 141 } 142 template <int N1_1, int N1_2, int N2_1, int N2_2> split_12mfem::StridedLayout2D143 static StridedLayout4D<N1_1,S1,N1_2,S1*N1_1,N2_1,S2,N2_2,S2*N2_1> split_12() 144 { 145 // S1*i1+S2*i2 == S1*(i1_1+N1_1*i1_2)+S2*(i2_1+N2_1*i2_2) 146 MFEM_STATIC_ASSERT(N1_1*N1_2 == N1 && N2_1*N2_2 == N2, 147 "invalid dimensions"); 148 return StridedLayout4D<N1_1,S1,N1_2,S1*N1_1,N2_1,S2,N2_2,S2*N2_1>(); 149 } merge_12mfem::StridedLayout2D150 static StridedLayout1D<N1*N2,(S1<S2)?S1:S2> merge_12() 151 { 152 // use: (S1*i1+S2*i2) == (S1*(i1+S2/S1*i2)) 153 // or (S1*i1+S2*i2) == (S2*(S1/S2*i1+i2)) 154 // assuming: S2 == S1*N1 || S1 == S2*N2 155 MFEM_STATIC_ASSERT(S2 == S1*N1 || S1 == S2*N2, "invalid reshape"); 156 return StridedLayout1D<N1*N2,(S1<S2)?S1:S2>(); 157 } transpose_12mfem::StridedLayout2D158 static StridedLayout2D<N2,S2,N1,S1> transpose_12() 159 { 160 return StridedLayout2D<N2,S2,N1,S1>(); 161 } 162 }; 163 164 template <int N1, int S1, int N2, int S2, int N3, int S3> 165 struct OffsetStridedLayout3D; 166 template <int N1, int S1, int N2, int S2, int N3, int S3, int N4, int S4> 167 struct OffsetStridedLayout4D; 168 169 template <int N1, int S1, int N2, int S2> 170 struct OffsetStridedLayout2D 171 { 172 static const int rank = 2; 173 static const int dim_1 = N1; 174 static const int dim_2 = N2; 175 static const int size = N1*N2; 176 177 int offset; 178 OffsetStridedLayout2Dmfem::OffsetStridedLayout2D179 OffsetStridedLayout2D() { } OffsetStridedLayout2Dmfem::OffsetStridedLayout2D180 OffsetStridedLayout2D(int offset_) : offset(offset_) { } indmfem::OffsetStridedLayout2D181 MFEM_HOST_DEVICE inline int ind(int i1, int i2) const 182 { 183 return offset+S1*i1+S2*i2; 184 } ind1mfem::OffsetStridedLayout2D185 OffsetStridedLayout1D<N2,S2> ind1(int i1) const 186 { 187 return OffsetStridedLayout1D<N2,S2>(offset+S1*i1); 188 } ind2mfem::OffsetStridedLayout2D189 OffsetStridedLayout1D<N1,S1> ind2(int i2) const 190 { 191 return OffsetStridedLayout1D<N1,S1>(offset+S2*i2); 192 } 193 194 template <int M1, int M2> submfem::OffsetStridedLayout2D195 OffsetStridedLayout2D<M1,S1,M2,S2> sub(int o1, int o2) const 196 { 197 return OffsetStridedLayout2D<M1,S1,M2,S2>(offset+S1*o1+S2*o2); 198 } 199 200 // reshape methods 201 202 template <int N1_1, int N1_2> split_1mfem::OffsetStridedLayout2D203 OffsetStridedLayout3D<N1_1,S1,N1_2,S1*N1_1,N2,S2> split_1() const 204 { 205 // S1*i1+S2*i2 == S1*(i1_1+N1_1*i1_2)+S2*i2 206 MFEM_STATIC_ASSERT(N1_1*N1_2 == N1, "invalid dimensions"); 207 return OffsetStridedLayout3D<N1_1,S1,N1_2,S1*N1_1,N2,S2>(offset); 208 } 209 template <int N2_1, int N2_2> split_2mfem::OffsetStridedLayout2D210 OffsetStridedLayout3D<N1,S1,N2_1,S2,N2_2,S2*N2_1> split_2() const 211 { 212 // S1*i1+S2*i2 == S1*i1+S2*(i2_1*N2_1*i2_2) 213 MFEM_STATIC_ASSERT(N2_1*N2_2 == N2, "invalid dimensions"); 214 return OffsetStridedLayout3D<N1,S1,N2_1,S2,N2_2,S2*N2_1>(offset); 215 } 216 template <int N1_1, int N1_2, int N2_1, int N2_2> 217 OffsetStridedLayout4D<N1_1,S1,N1_2,S1*N1_1,N2_1,S2,N2_2,S2*N2_1> split_12mfem::OffsetStridedLayout2D218 split_12() const 219 { 220 // S1*i1+S2*i2 == S1*(i1_1+N1_1*i1_2)+S2*(i2_1+N2_1*i2_2) 221 MFEM_STATIC_ASSERT(N1_1*N1_2 == N1 && N2_1*N2_2 == N2, 222 "invalid dimensions"); 223 return OffsetStridedLayout4D< 224 N1_1,S1,N1_2,S1*N1_1,N2_1,S2,N2_2,S2*N2_1>(offset); 225 } merge_12mfem::OffsetStridedLayout2D226 OffsetStridedLayout1D<N1*N2,(S1<S2)?S1:S2> merge_12() const 227 { 228 // use: (S1*i1+S2*i2) == (S1*(i1+S2/S1*i2)) 229 // or (S1*i1+S2*i2) == (S2*(S1/S2*i1+i2)) 230 // assuming: S2 == S1*N1 || S1 == S2*N2 231 MFEM_STATIC_ASSERT(S2 == S1*N1 || S1 == S2*N2, "invalid reshape"); 232 return OffsetStridedLayout1D<N1*N2,(S1<S2)?S1:S2>(offset); 233 } transpose_12mfem::OffsetStridedLayout2D234 OffsetStridedLayout2D<N2,S2,N1,S1> transpose_12() const 235 { 236 return OffsetStridedLayout2D<N2,S2,N1,S1>(offset); 237 } 238 }; 239 240 template <int N1, int S1, int N2, int S2, int N3, int S3> 241 struct StridedLayout3D 242 { 243 static const int rank = 3; 244 static const int dim_1 = N1; 245 static const int dim_2 = N2; 246 static const int dim_3 = N3; 247 static const int size = N1*N2*N3; 248 indmfem::StridedLayout3D249 static inline int ind(int i1, int i2, int i3) 250 { 251 return S1*i1+S2*i2+S3*i3; 252 } ind1mfem::StridedLayout3D253 static OffsetStridedLayout2D<N2,S2,N3,S3> ind1(int i1) 254 { 255 return OffsetStridedLayout2D<N2,S2,N3,S3>(S1*i1); 256 } ind2mfem::StridedLayout3D257 static OffsetStridedLayout2D<N1,S1,N3,S3> ind2(int i2) 258 { 259 return OffsetStridedLayout2D<N1,S1,N3,S3>(S2*i2); 260 } ind3mfem::StridedLayout3D261 static OffsetStridedLayout2D<N1,S1,N2,S2> ind3(int i3) 262 { 263 return OffsetStridedLayout2D<N1,S1,N2,S2>(S3*i3); 264 } 265 266 // reshape methods 267 merge_12mfem::StridedLayout3D268 static StridedLayout2D<N1*N2,S1,N3,S3> merge_12() 269 { 270 // use: (S1*i1+S2*i2+S3*i3) == (S1*(i1+S2/S1*i2)+S3*i3) 271 // assuming: S2 == S1*N1 272 MFEM_STATIC_ASSERT(S2 == S1*N1, "invalid reshape"); 273 return StridedLayout2D<N1*N2,S1,N3,S3>(); 274 // alternative: 275 // use: (S1*i1+S2*i2+S3*i3) == (S2*(S1/S2*i1+i2)+S3*i3) 276 // assuming: S1 == S2*N2 277 // result is: StridedLayout2D<N1*N2,S2,N3,S3> 278 } merge_23mfem::StridedLayout3D279 static StridedLayout2D<N1,S1,N2*N3,S2> merge_23() 280 { 281 // use: (S1*i1+S2*i2+S3*i3) == (S1*i1+S2*(i2+S3/S2*i3)) 282 // assuming: S3 == S2*N2 283 MFEM_STATIC_ASSERT(S3 == S2*N2, "invalid reshape"); 284 return StridedLayout2D<N1,S1,N2*N3,S2>(); 285 } 286 287 template <int N1_1, int N1_2> split_1mfem::StridedLayout3D288 static StridedLayout4D<N1_1,S1,N1_2,S1*N1_1,N2,S2,N3,S3> split_1() 289 { 290 // S1*i1+S2*i2+S3*i3 == S1*(i1_1+N1_1*i1_2)+S2*i2+S3*i3 291 MFEM_STATIC_ASSERT(N1_1*N1_2 == N1, "invalid dimensions"); 292 return StridedLayout4D<N1_1,S1,N1_2,S1*N1_1,N2,S2,N3,S3>(); 293 } 294 template <int N2_1, int N2_2> split_2mfem::StridedLayout3D295 static StridedLayout4D<N1,S1,N2_1,S2,N2_2,S2*N2_1,N3,S3> split_2() 296 { 297 // S1*i1+S2*i2+S3*i3 == S1*i1+S2*(i2_1+N2_1*i2_2)+S3*i3 298 MFEM_STATIC_ASSERT(N2_1*N2_2 == N2, "invalid dimensions"); 299 return StridedLayout4D<N1,S1,N2_1,S2,N2_2,S2*N2_1,N3,S3>(); 300 } 301 template <int N3_1, int N3_2> split_3mfem::StridedLayout3D302 static StridedLayout4D<N1,S1,N2,S2,N3_1,S3,N3_2,S3*N3_1> split_3() 303 { 304 // S1*i1+S2*i2+S3*i3 == S1*i1+S2*i2+S3*(i3_1+N3_1*i3_2) 305 MFEM_STATIC_ASSERT(N3_1*N3_2 == N3, "invalid dimensions"); 306 return StridedLayout4D<N1,S1,N2,S2,N3_1,S3,N3_2,S3*N3_1>(); 307 } 308 transpose_12mfem::StridedLayout3D309 static StridedLayout3D<N2,S2,N1,S1,N3,S3> transpose_12() 310 { 311 return StridedLayout3D<N2,S2,N1,S1,N3,S3>(); 312 } transpose_13mfem::StridedLayout3D313 static StridedLayout3D<N3,S3,N2,S2,N1,S1> transpose_13() 314 { 315 return StridedLayout3D<N3,S3,N2,S2,N1,S1>(); 316 } transpose_23mfem::StridedLayout3D317 static StridedLayout3D<N1,S1,N3,S3,N2,S2> transpose_23() 318 { 319 return StridedLayout3D<N1,S1,N3,S3,N2,S2>(); 320 } 321 }; 322 323 template <int N1, int S1, int N2, int S2, int N3, int S3> 324 struct OffsetStridedLayout3D 325 { 326 static const int rank = 3; 327 static const int dim_1 = N1; 328 static const int dim_2 = N2; 329 static const int dim_3 = N3; 330 static const int size = N1*N2*N3; 331 332 int offset; 333 OffsetStridedLayout3Dmfem::OffsetStridedLayout3D334 OffsetStridedLayout3D() { } OffsetStridedLayout3Dmfem::OffsetStridedLayout3D335 OffsetStridedLayout3D(int offset_) : offset(offset_) { } indmfem::OffsetStridedLayout3D336 inline int ind(int i1, int i2, int i3) const 337 { 338 return offset+S1*i1+S2*i2+S3*i3; 339 } ind1mfem::OffsetStridedLayout3D340 OffsetStridedLayout2D<N2,S2,N3,S3> ind1(int i1) const 341 { 342 return OffsetStridedLayout2D<N2,S2,N3,S3>(offset+S1*i1); 343 } ind2mfem::OffsetStridedLayout3D344 OffsetStridedLayout2D<N1,S1,N3,S3> ind2(int i2) const 345 { 346 return OffsetStridedLayout2D<N1,S1,N3,S3>(offset+S2*i2); 347 } ind3mfem::OffsetStridedLayout3D348 OffsetStridedLayout2D<N1,S1,N2,S2> ind3(int i3) const 349 { 350 return OffsetStridedLayout2D<N1,S1,N2,S2>(offset+S3*i3); 351 } 352 353 // reshape methods 354 merge_12mfem::OffsetStridedLayout3D355 OffsetStridedLayout2D<N1*N2,S1,N3,S3> merge_12() const 356 { 357 // use: (S1*i1+S2*i2+S3*i3) == (S1*(i1+S2/S1*i2)+S3*i3) 358 // assuming: S2 == S1*N1 359 MFEM_STATIC_ASSERT(S2 == S1*N1, "invalid reshape"); 360 return OffsetStridedLayout2D<N1*N2,S1,N3,S3>(offset); 361 } merge_23mfem::OffsetStridedLayout3D362 OffsetStridedLayout2D<N1,S1,N2*N3,S2> merge_23() const 363 { 364 // use: (S1*i1+S2*i2+S3*i3) == (S1*i1+S2*(i2+S3/S2*i3)) 365 // assuming: S3 == S2*N2 366 MFEM_STATIC_ASSERT(S3 == S2*N2, "invalid reshape"); 367 return OffsetStridedLayout2D<N1,S1,N2*N3,S2>(offset); 368 } 369 370 template <int N1_1, int N1_2> split_1mfem::OffsetStridedLayout3D371 OffsetStridedLayout4D<N1_1,S1,N1_2,S1*N1_1,N2,S2,N3,S3> split_1() const 372 { 373 // S1*i1+S2*i2+S3*i3 == S1*(i1_1+N1_1*i1_2)+S2*i2+S3*i3 374 MFEM_STATIC_ASSERT(N1_1*N1_2 == N1, "invalid dimensions"); 375 return OffsetStridedLayout4D<N1_1,S1,N1_2,S1*N1_1,N2,S2,N3,S3>(offset); 376 } 377 template <int N2_1, int N2_2> split_2mfem::OffsetStridedLayout3D378 OffsetStridedLayout4D<N1,S1,N2_1,S2,N2_2,S2*N2_1,N3,S3> split_2() const 379 { 380 // S1*i1+S2*i2+S3*i3 == S1*i1+S2*(i2_1+N2_1*i2_2)+S3*i3 381 MFEM_STATIC_ASSERT(N2_1*N2_2 == N2, "invalid dimensions"); 382 return OffsetStridedLayout4D<N1,S1,N2_1,S2,N2_2,S2*N2_1,N3,S3>(offset); 383 } 384 }; 385 386 template <int N1, int S1, int N2, int S2, int N3, int S3, int N4, int S4> 387 struct StridedLayout4D 388 { 389 static const int rank = 4; 390 static const int dim_1 = N1; 391 static const int dim_2 = N2; 392 static const int dim_3 = N3; 393 static const int dim_4 = N4; 394 static const int size = N1*N2*N3*N4; 395 indmfem::StridedLayout4D396 static inline int ind(int i1, int i2, int i3, int i4) 397 { 398 return S1*i1+S2*i2+S3*i3+S4*i4; 399 } ind23mfem::StridedLayout4D400 static OffsetStridedLayout2D<N1,S1,N4,S4> ind23(int i2, int i3) 401 { 402 return OffsetStridedLayout2D<N1,S1,N4,S4>(S2*i2+S3*i3); 403 } ind14mfem::StridedLayout4D404 static OffsetStridedLayout2D<N2,S2,N3,S3> ind14(int i1, int i4) 405 { 406 return OffsetStridedLayout2D<N2,S2,N3,S3>(S1*i1+S4*i4); 407 } ind4mfem::StridedLayout4D408 static OffsetStridedLayout3D<N1,S1,N2,S2,N3,S3> ind4(int i4) 409 { 410 return OffsetStridedLayout3D<N1,S1,N2,S2,N3,S3>(S4*i4); 411 } 412 merge_12mfem::StridedLayout4D413 static StridedLayout3D<N1*N2,S1,N3,S3,N4,S4> merge_12() 414 { 415 // use: (S1*i1+S2*i2+S3*i3+S4*i4) == (S1*(i1+S2/S1*i2)+S3*i3+S4*i4) 416 // assuming S2 == S1*N1 417 MFEM_STATIC_ASSERT(S2 == S1*N1, "invalid reshape"); 418 return StridedLayout3D<N1*N2,S1,N3,S3,N4,S4>(); 419 } merge_34mfem::StridedLayout4D420 static StridedLayout3D<N1,S1,N2,S2,N3*N4,S3> merge_34() 421 { 422 // use: (S1*i1+S2*i2+S3*i3+S4*i4) == (S1*i1+S2*i2+S3*(i3+S4/S3*i4)) 423 // assuming S4 == S3*N3 424 MFEM_STATIC_ASSERT(S4 == S3*N3, "invalid reshape"); 425 return StridedLayout3D<N1,S1,N2,S2,N3*N4,S3>(); 426 } 427 }; 428 429 template <int N1, int S1, int N2, int S2, int N3, int S3, int N4, int S4> 430 struct OffsetStridedLayout4D 431 { 432 static const int rank = 4; 433 static const int dim_1 = N1; 434 static const int dim_2 = N2; 435 static const int dim_3 = N3; 436 static const int dim_4 = N4; 437 static const int size = N1*N2*N3*N4; 438 439 int offset; 440 OffsetStridedLayout4Dmfem::OffsetStridedLayout4D441 OffsetStridedLayout4D() { } OffsetStridedLayout4Dmfem::OffsetStridedLayout4D442 OffsetStridedLayout4D(int offset_) : offset(offset_) { } indmfem::OffsetStridedLayout4D443 inline int ind(int i1, int i2, int i3, int i4) const 444 { 445 return offset+S1*i1+S2*i2+S3*i3+S4*i4; 446 } 447 }; 448 449 template <int N1, int N2> 450 struct ColumnMajorLayout2D 451 : public StridedLayout2D<N1,1,N2,N1> { }; 452 453 template <int N1, int N2, int N3> 454 struct ColumnMajorLayout3D 455 : public StridedLayout3D<N1,1,N2,N1,N3,N1*N2> { }; 456 457 template <int N1, int N2, int N3, int N4> 458 struct ColumnMajorLayout4D 459 : public StridedLayout4D<N1,1,N2,N1,N3,N1*N2,N4,N1*N2*N3> { }; 460 461 462 // Vector layout classes 463 464 class DynamicVectorLayout 465 { 466 public: 467 static const int vec_dim = 0; // 0 - dynamic 468 469 protected: 470 int scal_stride, comp_stride; 471 int num_components; 472 Init(Ordering::Type ordering,int scalar_size,int num_comp)473 void Init(Ordering::Type ordering, int scalar_size, int num_comp) 474 { 475 num_components = num_comp; 476 if (ordering == Ordering::byNODES) 477 { 478 scal_stride = 1; 479 comp_stride = scalar_size; 480 } 481 else 482 { 483 scal_stride = num_comp; 484 comp_stride = 1; 485 } 486 } 487 488 public: DynamicVectorLayout(Ordering::Type ordering,int scalar_size,int num_comp)489 DynamicVectorLayout(Ordering::Type ordering, int scalar_size, int num_comp) 490 { 491 Init(ordering, scalar_size, num_comp); 492 } DynamicVectorLayout(const FiniteElementSpace & fes)493 DynamicVectorLayout(const FiniteElementSpace &fes) 494 { 495 Init(fes.GetOrdering(), fes.GetNDofs(), fes.GetVDim()); 496 } 497 // default copy constructor 498 NumComponents() const499 int NumComponents() const { return num_components; } 500 ind(int scalar_idx,int comp_idx) const501 int ind(int scalar_idx, int comp_idx) const 502 { 503 return scal_stride * scalar_idx + comp_stride * comp_idx; 504 } 505 Matches(const FiniteElementSpace & fes)506 static bool Matches(const FiniteElementSpace &fes) 507 { 508 return true; 509 } 510 }; 511 512 // The default value (NumComp = 0) indicates that the number of components is 513 // dynamic, i.e. it will be specified at run-time. 514 template <Ordering::Type Ord, int NumComp = 0> 515 class VectorLayout 516 { 517 public: 518 static const int vec_dim = NumComp; 519 520 protected: 521 int num_components, scalar_size; 522 523 public: VectorLayout(int scalar_size_,int num_comp_=NumComp)524 VectorLayout(int scalar_size_, int num_comp_ = NumComp) 525 : num_components(num_comp_), 526 scalar_size(scalar_size_) 527 { 528 MFEM_ASSERT(NumComp == 0 || num_components == NumComp, 529 "invalid number of components"); 530 } 531 VectorLayout(const FiniteElementSpace & fes)532 VectorLayout(const FiniteElementSpace &fes) 533 : num_components(fes.GetVDim()), 534 scalar_size(fes.GetNDofs()) 535 { 536 MFEM_ASSERT(fes.GetOrdering() == Ord, "ordering mismatch"); 537 MFEM_ASSERT(NumComp == 0 || num_components == NumComp, 538 "invalid number of components"); 539 } 540 // default copy constructor 541 NumComponents() const542 int NumComponents() const { return (NumComp ? NumComp : num_components); } 543 ind(int scalar_idx,int comp_idx) const544 int ind(int scalar_idx, int comp_idx) const 545 { 546 if (Ord == Ordering::byNODES) 547 { 548 return scalar_idx + comp_idx * scalar_size; 549 } 550 else 551 { 552 return comp_idx + (NumComp ? NumComp : num_components) * scalar_idx; 553 } 554 } 555 Matches(const FiniteElementSpace & fes)556 static bool Matches(const FiniteElementSpace &fes) 557 { 558 return (Ord == fes.GetOrdering() && 559 (NumComp == 0 || NumComp == fes.GetVDim())); 560 } 561 }; 562 563 class ScalarLayout 564 { 565 public: 566 static const int vec_dim = 1; 567 ScalarLayout()568 ScalarLayout() { } 569 ScalarLayout(const FiniteElementSpace & fes)570 ScalarLayout(const FiniteElementSpace &fes) 571 { 572 MFEM_ASSERT(fes.GetVDim() == 1, "invalid number of components"); 573 } 574 NumComponents() const575 int NumComponents() const { return 1; } 576 ind(int scalar_idx,int comp_idx) const577 int ind(int scalar_idx, int comp_idx) const { return scalar_idx; } 578 Matches(const FiniteElementSpace & fes)579 static bool Matches(const FiniteElementSpace &fes) 580 { 581 return (fes.GetVDim() == 1); 582 } 583 }; 584 585 } // namespace mfem 586 587 #endif // MFEM_TEMPLATE_LAYOUT 588