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_NURBS
13 #define MFEM_NURBS
14
15 #include "../config/config.hpp"
16 #include "../general/table.hpp"
17 #include "../linalg/vector.hpp"
18 #include "element.hpp"
19 #include "mesh.hpp"
20 #ifdef MFEM_USE_MPI
21 #include "../general/communication.hpp"
22 #endif
23 #include <iostream>
24
25 namespace mfem
26 {
27
28 class GridFunction;
29
30
31 class KnotVector
32 {
33 protected:
34 static const int MaxOrder;
35
36 Vector knot;
37 int Order, NumOfControlPoints, NumOfElements;
38
39 public:
40 /// Create KnotVector
KnotVector()41 KnotVector() { }
42 KnotVector(std::istream &input);
43 KnotVector(int Order_, int NCP);
KnotVector(const KnotVector & kv)44 KnotVector(const KnotVector &kv) { (*this) = kv; }
45
46 KnotVector &operator=(const KnotVector &kv);
47
GetNE() const48 int GetNE() const { return NumOfElements; }
GetNKS() const49 int GetNKS() const { return NumOfControlPoints - Order; }
GetNCP() const50 int GetNCP() const { return NumOfControlPoints; }
GetOrder() const51 int GetOrder() const { return Order; }
Size() const52 int Size() const { return knot.Size(); }
53
54 /// Count the number of elements
55 void GetElements();
56
isElement(int i) const57 bool isElement(int i) const { return (knot(Order+i) != knot(Order+i+1)); }
58
getKnotLocation(double xi,int ni) const59 double getKnotLocation(double xi, int ni) const
60 { return (xi*knot(ni+1) + (1. - xi)*knot(ni)); }
61
62 int findKnotSpan(double u) const;
63
64 void CalcShape (Vector &shape, int i, double xi) const;
65 void CalcDShape (Vector &grad, int i, double xi) const;
66 void CalcDnShape(Vector &gradn, int n, int i, double xi) const;
CalcD2Shape(Vector & grad2,int i,double xi) const67 void CalcD2Shape(Vector &grad2, int i, double xi) const
68 { CalcDnShape(grad2, 2, i, xi); }
69
70 void Difference(const KnotVector &kv, Vector &diff) const;
71 void UniformRefinement(Vector &newknots) const;
72 /** Return a new KnotVector with elevated degree by repeating the endpoints
73 of the knot vector. */
74 KnotVector *DegreeElevate(int t) const;
75
76 void Flip();
77
78 void Print(std::ostream &out) const;
79
80 void PrintFunctions(std::ostream &out, int samples=11) const;
81
82 /// Destroys KnotVector
~KnotVector()83 ~KnotVector() { }
84
operator [](int i)85 double &operator[](int i) { return knot(i); }
operator [](int i) const86 const double &operator[](int i) const { return knot(i); }
87 };
88
89
90 class NURBSPatch
91 {
92 protected:
93 int ni, nj, nk, Dim;
94 double *data; // the layout of data is: (Dim x ni x nj x nk)
95
96 Array<KnotVector *> kv;
97
98 int sd, nd;
99
100 void swap(NURBSPatch *np);
101
102 // Special B-NET access functions
103 int SetLoopDirection(int dir);
104 inline double &operator()(int i, int j);
105 inline const double &operator()(int i, int j) const;
106
107 void init(int dim_);
108
109 NURBSPatch(NURBSPatch *parent, int dir, int Order, int NCP);
110
111 public:
112 NURBSPatch(const NURBSPatch &orig);
113 NURBSPatch(std::istream &input);
114 NURBSPatch(const KnotVector *kv0, const KnotVector *kv1, int dim_);
115 NURBSPatch(const KnotVector *kv0, const KnotVector *kv1,
116 const KnotVector *kv2, int dim_);
117 NURBSPatch(Array<const KnotVector *> &kv, int dim_);
118
119 ~NURBSPatch();
120
121 void Print(std::ostream &out) const;
122
123 void DegreeElevate(int dir, int t);
124 void KnotInsert (int dir, const KnotVector &knot);
125 void KnotInsert (int dir, const Vector &knot);
126
127 void KnotInsert(Array<Vector *> &knot);
128 void KnotInsert(Array<KnotVector *> &knot);
129
130 void DegreeElevate(int t);
131 void UniformRefinement();
132
133 // Return the number of components stored in the NURBSPatch
GetNC() const134 int GetNC() const { return Dim; }
GetNKV() const135 int GetNKV() const { return kv.Size(); }
GetKV(int i)136 KnotVector *GetKV(int i) { return kv[i]; }
137
138 // Standard B-NET access functions
139 inline double &operator()(int i, int j, int l);
140 inline const double &operator()(int i, int j, int l) const;
141
142 inline double &operator()(int i, int j, int k, int l);
143 inline const double &operator()(int i, int j, int k, int l) const;
144
145 static void Get3DRotationMatrix(double n[], double angle, double r,
146 DenseMatrix &T);
147 void FlipDirection(int dir);
148 void SwapDirections(int dir1, int dir2);
149 void Rotate3D(double normal[], double angle);
150 int MakeUniformDegree(int degree = -1);
151 friend NURBSPatch *Interpolate(NURBSPatch &p1, NURBSPatch &p2);
152 friend NURBSPatch *Revolve3D(NURBSPatch &patch, double n[], double ang,
153 int times);
154 };
155
156
157 #ifdef MFEM_USE_MPI
158 class ParNURBSExtension;
159 #endif
160
161 class NURBSPatchMap;
162
163
164 class NURBSExtension
165 {
166 #ifdef MFEM_USE_MPI
167 friend class ParNURBSExtension;
168 #endif
169 friend class NURBSPatchMap;
170
171 protected:
172 int mOrder; // see GetOrder() for description
173 Array<int> mOrders;
174 int NumOfKnotVectors;
175 // global entity counts
176 int NumOfVertices, NumOfElements, NumOfBdrElements, NumOfDofs;
177 // local entity counts
178 int NumOfActiveVertices, NumOfActiveElems, NumOfActiveBdrElems;
179 int NumOfActiveDofs;
180
181 Array<int> activeVert; // activeVert[glob_vert] = loc_vert or -1
182 Array<bool> activeElem;
183 Array<bool> activeBdrElem;
184 Array<int> activeDof; // activeDof[glob_dof] = loc_dof + 1 or 0
185
186 Mesh *patchTopo;
187 int own_topo;
188 Array<int> edge_to_knot;
189 Array<KnotVector *> knotVectors;
190 Vector weights;
191
192 // periodic BC info:
193 // - dof 2 dof map
194 // - master and slave boundary indices
195 Array<int> d_to_d;
196 Array<int> master;
197 Array<int> slave;
198
199 // global offsets, meshOffsets == meshVertexOffsets
200 Array<int> v_meshOffsets;
201 Array<int> e_meshOffsets;
202 Array<int> f_meshOffsets;
203 Array<int> p_meshOffsets;
204
205 // global offsets, spaceOffsets == dofOffsets
206 Array<int> v_spaceOffsets;
207 Array<int> e_spaceOffsets;
208 Array<int> f_spaceOffsets;
209 Array<int> p_spaceOffsets;
210
211 Table *el_dof, *bel_dof;
212
213 Array<int> el_to_patch;
214 Array<int> bel_to_patch;
215 Array2D<int> el_to_IJK; // IJK are "knot-span" indices!
216 Array2D<int> bel_to_IJK; // they are NOT element indices!
217
218 Array<NURBSPatch *> patches;
219
220 inline int KnotInd(int edge) const;
221 inline KnotVector *KnotVec(int edge);
222 inline const KnotVector *KnotVec(int edge) const;
223 inline const KnotVector *KnotVec(int edge, int oedge, int *okv) const;
224
225 void CheckPatches();
226 void CheckBdrPatches();
227
228 void GetPatchKnotVectors (int p, Array<KnotVector *> &kv);
229 void GetPatchKnotVectors (int p, Array<const KnotVector *> &kv) const;
230 void GetBdrPatchKnotVectors(int p, Array<KnotVector *> &kv);
231 void GetBdrPatchKnotVectors(int p, Array<const KnotVector *> &kv) const;
232
233 void SetOrderFromOrders();
234 void SetOrdersFromKnotVectors();
235
236 // periodic BC helper functions
237 void InitDofMap();
238 void ConnectBoundaries();
239 void ConnectBoundaries2D(int bnd0, int bnd1);
240 void ConnectBoundaries3D(int bnd0, int bnd1);
DofMap(int dof) const241 int DofMap(int dof) const
242 {
243 return (d_to_d.Size() > 0 )? d_to_d[dof] : dof;
244 };
245
246 // also count the global NumOfVertices and the global NumOfDofs
247 void GenerateOffsets();
248 // count the global NumOfElements
249 void CountElements();
250 // count the global NumOfBdrElements
251 void CountBdrElements();
252
253 // generate the mesh elements
254 void Get2DElementTopo(Array<Element *> &elements) const;
255 void Get3DElementTopo(Array<Element *> &elements) const;
256
257 // generate the boundary mesh elements
258 void Get2DBdrElementTopo(Array<Element *> &boundary) const;
259 void Get3DBdrElementTopo(Array<Element *> &boundary) const;
260
261
262 // FE space generation functions
263
264 // based on activeElem, count NumOfActiveDofs, generate el_dof,
265 // el_to_patch, el_to_IJK, activeDof map (global-to-local)
266 void GenerateElementDofTable();
267
268 // generate elem_to_global-dof table for the active elements
269 // define el_to_patch, el_to_IJK, activeDof (as bool)
270 void Generate2DElementDofTable();
271 void Generate3DElementDofTable();
272
273 // call after GenerateElementDofTable
274 void GenerateBdrElementDofTable();
275
276 // generate the bdr-elem_to_global-dof table for the active bdr. elements
277 // define bel_to_patch, bel_to_IJK
278 void Generate2DBdrElementDofTable();
279 void Generate3DBdrElementDofTable();
280
281 // FE --> Patch translation functions
282 void GetPatchNets (const Vector &Nodes, int vdim);
283 void Get2DPatchNets(const Vector &Nodes, int vdim);
284 void Get3DPatchNets(const Vector &Nodes, int vdim);
285
286 // Patch --> FE translation functions
287 // Side effects: delete the patches, update the weights from the patches
288 void SetSolutionVector (Vector &Nodes, int vdim);
289 void Set2DSolutionVector(Vector &Nodes, int vdim);
290 void Set3DSolutionVector(Vector &Nodes, int vdim);
291
292 // determine activeVert, NumOfActiveVertices from the activeElem array
293 void GenerateActiveVertices();
294
295 // determine activeBdrElem, NumOfActiveBdrElems
296 void GenerateActiveBdrElems();
297
298 void MergeWeights(Mesh *mesh_array[], int num_pieces);
299
300 // to be used by ParNURBSExtension constructor(s)
NURBSExtension()301 NURBSExtension() { }
302
303 public:
304 /// Copy constructor: deep copy
305 NURBSExtension(const NURBSExtension &orig);
306 /// Read-in a NURBSExtension
307 NURBSExtension(std::istream &input);
308 /** @brief Create a NURBSExtension with elevated order by repeating the
309 endpoints of the knot vectors and using uniform weights of 1. */
310 /** If a knot vector in @a parent already has order greater than or equal to
311 @a newOrder, it will be used unmodified. */
312 NURBSExtension(NURBSExtension *parent, int newOrder);
313 /** @brief Create a NURBSExtension with elevated knot vector orders (by
314 repeating the endpoints of the knot vectors and using uniform weights of
315 1) as given by the array @a newOrders. */
316 /** If a knot vector in @a parent already has order greater than or equal to
317 the corresponding entry in @a newOrder, it will be used unmodified. */
318 NURBSExtension(NURBSExtension *parent, const Array<int> &newOrders);
319 /// Construct a NURBSExtension by merging a partitioned NURBS mesh
320 NURBSExtension(Mesh *mesh_array[], int num_pieces);
321
322 // Generate connections between boundaries, such as periodic BCs
323 void ConnectBoundaries(Array<int> &master, Array<int> &slave);
GetMaster() const324 const Array<int> &GetMaster() const { return master; };
GetMaster()325 Array<int> &GetMaster() { return master; };
GetSlave() const326 const Array<int> &GetSlave() const { return slave; };
GetSlave()327 Array<int> &GetSlave() { return slave; };
328 void MergeGridFunctions(GridFunction *gf_array[], int num_pieces,
329 GridFunction &merged);
330
331 /// Destroy a NURBSExtension
332 virtual ~NURBSExtension();
333
334 // Print functions
335 void Print(std::ostream &out) const;
336 void PrintCharacteristics(std::ostream &out) const;
337 void PrintFunctions(const char *filename, int samples=11) const;
338
339 // Meta data functions
Dimension() const340 int Dimension() const { return patchTopo->Dimension(); }
GetNP() const341 int GetNP() const { return patchTopo->GetNE(); }
GetNBP() const342 int GetNBP() const { return patchTopo->GetNBE(); }
343
344 /// Read-only access to the orders of all knot vectors.
GetOrders() const345 const Array<int> &GetOrders() const { return mOrders; }
346 /** @brief If all orders are identical, return that number. Otherwise, return
347 NURBSFECollection::VariableOrder. */
GetOrder() const348 int GetOrder() const { return mOrder; }
349
GetNKV() const350 int GetNKV() const { return NumOfKnotVectors; }
351
GetGNV() const352 int GetGNV() const { return NumOfVertices; }
GetNV() const353 int GetNV() const { return NumOfActiveVertices; }
GetGNE() const354 int GetGNE() const { return NumOfElements; }
GetNE() const355 int GetNE() const { return NumOfActiveElems; }
GetGNBE() const356 int GetGNBE() const { return NumOfBdrElements; }
GetNBE() const357 int GetNBE() const { return NumOfActiveBdrElems; }
358
GetNTotalDof() const359 int GetNTotalDof() const { return NumOfDofs; }
GetNDof() const360 int GetNDof() const { return NumOfActiveDofs; }
361
362 // Knotvector read-only access function
GetKnotVector(int i) const363 const KnotVector *GetKnotVector(int i) const { return knotVectors[i]; }
364
365 // Mesh generation functions
366 void GetElementTopo (Array<Element *> &elements) const;
367 void GetBdrElementTopo(Array<Element *> &boundary) const;
368
HavePatches() const369 bool HavePatches() const { return (patches.Size() != 0); }
370
GetElementDofTable()371 Table *GetElementDofTable() { return el_dof; }
GetBdrElementDofTable()372 Table *GetBdrElementDofTable() { return bel_dof; }
373
374 void GetVertexLocalToGlobal(Array<int> &lvert_vert);
375 void GetElementLocalToGlobal(Array<int> &lelem_elem);
376
377 // Load functions
378 void LoadFE(int i, const FiniteElement *FE) const;
379 void LoadBE(int i, const FiniteElement *BE) const;
380
GetWeights() const381 const Vector &GetWeights() const { return weights; }
GetWeights()382 Vector &GetWeights() { return weights; }
383
384 // Translation functions: from FE coordinates to IJK patch
385 // format and vice versa
386 void ConvertToPatches(const Vector &Nodes);
387 void SetKnotsFromPatches();
388 void SetCoordsFromPatches(Vector &Nodes);
389
390 // Read a GridFunction written patch-by-patch, e.g. with PrintSolution().
391 void LoadSolution(std::istream &input, GridFunction &sol) const;
392 // Write a GridFunction patch-by-patch.
393 void PrintSolution(const GridFunction &sol, std::ostream &out) const;
394
395 // Refinement methods
396 // new_degree = max(old_degree, min(old_degree + rel_degree, degree))
397 void DegreeElevate(int rel_degree, int degree = 16);
398 void UniformRefinement();
399 void KnotInsert(Array<KnotVector *> &kv);
400 void KnotInsert(Array<Vector *> &kv);
401 };
402
403
404 #ifdef MFEM_USE_MPI
405 class ParNURBSExtension : public NURBSExtension
406 {
407 private:
408 int *partitioning;
409
410 Table *GetGlobalElementDofTable();
411 Table *Get2DGlobalElementDofTable();
412 Table *Get3DGlobalElementDofTable();
413
414 void SetActive(const int *partitioning, const Array<bool> &active_bel);
415 void BuildGroups(const int *partitioning, const Table &elem_dof);
416
417 public:
418 GroupTopology gtopo;
419
420 Array<int> ldof_group;
421
422 ParNURBSExtension(const ParNURBSExtension &orig);
423
424 ParNURBSExtension(MPI_Comm comm, NURBSExtension *parent, int *partitioning,
425 const Array<bool> &active_bel);
426
427 // Create a parallel version of 'parent' with partitioning as in
428 // 'par_parent'; the 'parent' object is destroyed.
429 // The 'parent' can be either a local NURBSExtension or a global one.
430 ParNURBSExtension(NURBSExtension *parent,
431 const ParNURBSExtension *par_parent);
432
~ParNURBSExtension()433 virtual ~ParNURBSExtension() { delete [] partitioning; }
434 };
435 #endif
436
437
438 class NURBSPatchMap
439 {
440 private:
441 const NURBSExtension *Ext;
442
443 int I, J, K, pOffset, opatch;
444 Array<int> verts, edges, faces, oedge, oface;
445
F(const int n,const int N)446 inline static int F(const int n, const int N)
447 { return (n < 0) ? 0 : ((n >= N) ? 2 : 1); }
448
Or1D(const int n,const int N,const int Or)449 inline static int Or1D(const int n, const int N, const int Or)
450 { return (Or > 0) ? n : (N - 1 - n); }
451
452 inline static int Or2D(const int n1, const int n2,
453 const int N1, const int N2, const int Or);
454
455 // also set verts, edges, faces, orientations etc
456 void GetPatchKnotVectors (int p, const KnotVector *kv[]);
457 void GetBdrPatchKnotVectors(int p, const KnotVector *kv[], int *okv);
458
459 public:
NURBSPatchMap(const NURBSExtension * ext)460 NURBSPatchMap(const NURBSExtension *ext) { Ext = ext; }
461
nx()462 int nx() { return I + 1; }
ny()463 int ny() { return J + 1; }
nz()464 int nz() { return K + 1; }
465
466 void SetPatchVertexMap(int p, const KnotVector *kv[]);
467 void SetPatchDofMap (int p, const KnotVector *kv[]);
468
469 void SetBdrPatchVertexMap(int p, const KnotVector *kv[], int *okv);
470 void SetBdrPatchDofMap (int p, const KnotVector *kv[], int *okv);
471
472 inline int operator()(const int i) const;
operator [](const int i) const473 inline int operator[](const int i) const { return (*this)(i); }
474
475 inline int operator()(const int i, const int j) const;
476
477 inline int operator()(const int i, const int j, const int k) const;
478 };
479
480
481 // Inline function implementations
482
operator ()(int i,int j)483 inline double &NURBSPatch::operator()(int i, int j)
484 {
485 return data[j%sd + sd*(i + (j/sd)*nd)];
486 }
487
operator ()(int i,int j) const488 inline const double &NURBSPatch::operator()(int i, int j) const
489 {
490 return data[j%sd + sd*(i + (j/sd)*nd)];
491 }
492
operator ()(int i,int j,int l)493 inline double &NURBSPatch::operator()(int i, int j, int l)
494 {
495 #ifdef MFEM_DEBUG
496 if (data == 0 || i < 0 || i >= ni || j < 0 || j >= nj || nk > 0 ||
497 l < 0 || l >= Dim)
498 {
499 mfem_error("NURBSPatch::operator() 2D");
500 }
501 #endif
502
503 return data[(i+j*ni)*Dim+l];
504 }
505
operator ()(int i,int j,int l) const506 inline const double &NURBSPatch::operator()(int i, int j, int l) const
507 {
508 #ifdef MFEM_DEBUG
509 if (data == 0 || i < 0 || i >= ni || j < 0 || j >= nj || nk > 0 ||
510 l < 0 || l >= Dim)
511 {
512 mfem_error("NURBSPatch::operator() const 2D");
513 }
514 #endif
515
516 return data[(i+j*ni)*Dim+l];
517 }
518
operator ()(int i,int j,int k,int l)519 inline double &NURBSPatch::operator()(int i, int j, int k, int l)
520 {
521 #ifdef MFEM_DEBUG
522 if (data == 0 || i < 0 || i >= ni || j < 0 || j >= nj || k < 0 ||
523 k >= nk || l < 0 || l >= Dim)
524 {
525 mfem_error("NURBSPatch::operator() 3D");
526 }
527 #endif
528
529 return data[(i+(j+k*nj)*ni)*Dim+l];
530 }
531
operator ()(int i,int j,int k,int l) const532 inline const double &NURBSPatch::operator()(int i, int j, int k, int l) const
533 {
534 #ifdef MFEM_DEBUG
535 if (data == 0 || i < 0 || i >= ni || j < 0 || j >= nj || k < 0 ||
536 k >= nk || l < 0 || l >= Dim)
537 {
538 mfem_error("NURBSPatch::operator() const 3D");
539 }
540 #endif
541
542 return data[(i+(j+k*nj)*ni)*Dim+l];
543 }
544
545
KnotInd(int edge) const546 inline int NURBSExtension::KnotInd(int edge) const
547 {
548 int kv = edge_to_knot[edge];
549 return (kv >= 0) ? kv : (-1-kv);
550 }
551
KnotVec(int edge)552 inline KnotVector *NURBSExtension::KnotVec(int edge)
553 {
554 return knotVectors[KnotInd(edge)];
555 }
556
KnotVec(int edge) const557 inline const KnotVector *NURBSExtension::KnotVec(int edge) const
558 {
559 return knotVectors[KnotInd(edge)];
560 }
561
KnotVec(int edge,int oedge,int * okv) const562 inline const KnotVector *NURBSExtension::KnotVec(int edge, int oedge, int *okv)
563 const
564 {
565 int kv = edge_to_knot[edge];
566 if (kv >= 0)
567 {
568 *okv = oedge;
569 return knotVectors[kv];
570 }
571 else
572 {
573 *okv = -oedge;
574 return knotVectors[-1-kv];
575 }
576 }
577
578
579 // static method
Or2D(const int n1,const int n2,const int N1,const int N2,const int Or)580 inline int NURBSPatchMap::Or2D(const int n1, const int n2,
581 const int N1, const int N2, const int Or)
582 {
583 // Needs testing
584 switch (Or)
585 {
586 case 0: return n1 + n2*N1;
587 case 1: return n2 + n1*N2;
588 case 2: return n2 + (N1 - 1 - n1)*N2;
589 case 3: return (N1 - 1 - n1) + n2*N1;
590 case 4: return (N1 - 1 - n1) + (N2 - 1 - n2)*N1;
591 case 5: return (N2 - 1 - n2) + (N1 - 1 - n1)*N2;
592 case 6: return (N2 - 1 - n2) + n1*N2;
593 case 7: return n1 + (N2 - 1 - n2)*N1;
594 }
595 #ifdef MFEM_DEBUG
596 mfem_error("NURBSPatchMap::Or2D");
597 #endif
598 return -1;
599 }
600
operator ()(const int i) const601 inline int NURBSPatchMap::operator()(const int i) const
602 {
603 const int i1 = i - 1;
604 switch (F(i1, I))
605 {
606 case 0: return verts[0];
607 case 1: return pOffset + Or1D(i1, I, opatch);
608 case 2: return verts[1];
609 }
610 #ifdef MFEM_DEBUG
611 mfem_error("NURBSPatchMap::operator() const 1D");
612 #endif
613 return -1;
614 }
615
operator ()(const int i,const int j) const616 inline int NURBSPatchMap::operator()(const int i, const int j) const
617 {
618 const int i1 = i - 1, j1 = j - 1;
619 switch (3*F(j1, J) + F(i1, I))
620 {
621 case 0: return verts[0];
622 case 1: return edges[0] + Or1D(i1, I, oedge[0]);
623 case 2: return verts[1];
624 case 3: return edges[3] + Or1D(j1, J, -oedge[3]);
625 case 4: return pOffset + Or2D(i1, j1, I, J, opatch);
626 case 5: return edges[1] + Or1D(j1, J, oedge[1]);
627 case 6: return verts[3];
628 case 7: return edges[2] + Or1D(i1, I, -oedge[2]);
629 case 8: return verts[2];
630 }
631 #ifdef MFEM_DEBUG
632 mfem_error("NURBSPatchMap::operator() const 2D");
633 #endif
634 return -1;
635 }
636
operator ()(const int i,const int j,const int k) const637 inline int NURBSPatchMap::operator()(const int i, const int j, const int k)
638 const
639 {
640 // Needs testing
641 const int i1 = i - 1, j1 = j - 1, k1 = k - 1;
642 switch (3*(3*F(k1, K) + F(j1, J)) + F(i1, I))
643 {
644 case 0: return verts[0];
645 case 1: return edges[0] + Or1D(i1, I, oedge[0]);
646 case 2: return verts[1];
647 case 3: return edges[3] + Or1D(j1, J, oedge[3]);
648 case 4: return faces[0] + Or2D(i1, J - 1 - j1, I, J, oface[0]);
649 case 5: return edges[1] + Or1D(j1, J, oedge[1]);
650 case 6: return verts[3];
651 case 7: return edges[2] + Or1D(i1, I, oedge[2]);
652 case 8: return verts[2];
653 case 9: return edges[8] + Or1D(k1, K, oedge[8]);
654 case 10: return faces[1] + Or2D(i1, k1, I, K, oface[1]);
655 case 11: return edges[9] + Or1D(k1, K, oedge[9]);
656 case 12: return faces[4] + Or2D(J - 1 - j1, k1, J, K, oface[4]);
657 case 13: return pOffset + I*(J*k1 + j1) + i1;
658 case 14: return faces[2] + Or2D(j1, k1, J, K, oface[2]);
659 case 15: return edges[11] + Or1D(k1, K, oedge[11]);
660 case 16: return faces[3] + Or2D(I - 1 - i1, k1, I, K, oface[3]);
661 case 17: return edges[10] + Or1D(k1, K, oedge[10]);
662 case 18: return verts[4];
663 case 19: return edges[4] + Or1D(i1, I, oedge[4]);
664 case 20: return verts[5];
665 case 21: return edges[7] + Or1D(j1, J, oedge[7]);
666 case 22: return faces[5] + Or2D(i1, j1, I, J, oface[5]);
667 case 23: return edges[5] + Or1D(j1, J, oedge[5]);
668 case 24: return verts[7];
669 case 25: return edges[6] + Or1D(i1, I, oedge[6]);
670 case 26: return verts[6];
671 }
672 #ifdef MFEM_DEBUG
673 mfem_error("NURBSPatchMap::operator() const 3D");
674 #endif
675 return -1;
676 }
677
678 }
679
680 #endif
681