MFEM  v4.6.0
Finite element discretization library
nurbs.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 #include <set>
25 
26 namespace mfem
27 {
28 
29 class GridFunction;
30 
31 
33 {
34 protected:
35  static const int MaxOrder;
36 
39 
40 public:
41  /// Create KnotVector
42  KnotVector() { }
43  KnotVector(std::istream &input);
44  KnotVector(int Order_, int NCP);
45  KnotVector(const KnotVector &kv) { (*this) = kv; }
46 
47  KnotVector &operator=(const KnotVector &kv);
48 
49  int GetNE() const { return NumOfElements; }
50  int GetNKS() const { return NumOfControlPoints - Order; }
51  int GetNCP() const { return NumOfControlPoints; }
52  int GetOrder() const { return Order; }
53  int Size() const { return knot.Size(); }
54 
55  /// Count the number of elements
56  void GetElements();
57 
58  bool isElement(int i) const { return (knot(Order+i) != knot(Order+i+1)); }
59 
60  double getKnotLocation(double xi, int ni) const
61  { return (xi*knot(ni+1) + (1. - xi)*knot(ni)); }
62 
63  int findKnotSpan(double u) const;
64 
65  void CalcShape (Vector &shape, int i, double xi) const;
66  void CalcDShape (Vector &grad, int i, double xi) const;
67  void CalcDnShape(Vector &gradn, int n, int i, double xi) const;
68  void CalcD2Shape(Vector &grad2, int i, double xi) const
69  { CalcDnShape(grad2, 2, i, xi); }
70 
71  /** Gives the locations of the maxima of the knotvector in reference space. The
72  function gives the knotspan @a ks, the coordinate in the knotspan @a xi
73  and the coordinate of the maximum in parameter space @a u */
74  void FindMaxima(Array<int> &ks, Vector &xi, Vector &u);
75  /** Global curve interpolation through the points @a x. @a x is an array with the
76  length of the spatial dimension containing vectors with spatial coordinates. The
77  controlpoints of the interpolated curve are given in @a x in the same form.*/
79 
80  void Difference(const KnotVector &kv, Vector &diff) const;
81  void UniformRefinement(Vector &newknots) const;
82  /** Return a new KnotVector with elevated degree by repeating the endpoints
83  of the knot vector. */
84  /// @note The returned object should be deleted by the caller.
85  KnotVector *DegreeElevate(int t) const;
86 
87  void Flip();
88 
89  void Print(std::ostream &out) const;
90 
91  void PrintFunctions(std::ostream &out, int samples=11) const;
92 
93  /// Destroys KnotVector
95 
96  double &operator[](int i) { return knot(i); }
97  const double &operator[](int i) const { return knot(i); }
98 };
99 
100 
102 {
103 protected:
104  int ni, nj, nk, Dim;
105  double *data; // the layout of data is: (Dim x ni x nj x nk)
106 
108 
109  // Special B-NET access functions
110  // - SetLoopDirection(int dir) flattens the multi-dimensional B-NET in the
111  // requested direction. It effectively creates a 1D net.
112  // - The slice(int, int) operator is the access function in that flattened structure.
113  // The first int gives the slice and the second int the element in that slice.
114  // - Both routines are used in 'InsertKnot', 'DegreeElevate' and 'UniformRefinement'.
115  // - In older implementations slice(int int) was implemented as operator()(int, int)
116  int nd; // Number of knots in flattened structure
117  int ls; // Number of variables per knot in flattened structure
118  int sd; // Stride for data access
119  int SetLoopDirection(int dir);
120  inline double &slice(int i, int j);
121  inline const double &slice(int i, int j) const;
122 
123  NURBSPatch(NURBSPatch *parent, int dir, int Order, int NCP);
124  void swap(NURBSPatch *np);
125  void init(int dim_);
126 
127 public:
128  NURBSPatch(const NURBSPatch &orig);
129  NURBSPatch(std::istream &input);
130  NURBSPatch(const KnotVector *kv0, const KnotVector *kv1, int dim_);
131  NURBSPatch(const KnotVector *kv0, const KnotVector *kv1,
132  const KnotVector *kv2, int dim_);
134 
135  NURBSPatch& operator=(const NURBSPatch&) = delete;
136 
137  ~NURBSPatch();
138 
139  void Print(std::ostream &out) const;
140 
141  void DegreeElevate(int dir, int t);
142  void KnotInsert (int dir, const KnotVector &knot);
143  void KnotInsert (int dir, const Vector &knot);
144 
145  void KnotInsert(Array<Vector *> &knot);
146  void KnotInsert(Array<KnotVector *> &knot);
147 
148  void DegreeElevate(int t);
149  void UniformRefinement();
150 
151  // Return the number of components stored in the NURBSPatch
152  int GetNC() const { return Dim; }
153  int GetNKV() const { return kv.Size(); }
154  /// @note The returned object should NOT be deleted by the caller.
155  KnotVector *GetKV(int i) { return kv[i]; }
156 
157  // Standard B-NET access functions
158  inline double &operator()(int i, int j);
159  inline const double &operator()(int i, int j) const;
160 
161  inline double &operator()(int i, int j, int l);
162  inline const double &operator()(int i, int j, int l) const;
163 
164  inline double &operator()(int i, int j, int k, int l);
165  inline const double &operator()(int i, int j, int k, int l) const;
166 
167  static void Get2DRotationMatrix(double angle,
168  DenseMatrix &T);
169  static void Get3DRotationMatrix(double n[], double angle, double r,
170  DenseMatrix &T);
171  void FlipDirection(int dir);
172  void SwapDirections(int dir1, int dir2);
173 
174  /// Rotate the NURBSPatch.
175  /** A rotation of a 2D NURBS-patch requires an angle only. Rotating
176  a 3D NURBS-patch requires a normal as well.*/
177  void Rotate(double angle, double normal[]= NULL);
178  void Rotate2D(double angle);
179  void Rotate3D(double normal[], double angle);
180 
181  int MakeUniformDegree(int degree = -1);
182  /// @note The returned object should be deleted by the caller.
183  friend NURBSPatch *Interpolate(NURBSPatch &p1, NURBSPatch &p2);
184  /// @note The returned object should be deleted by the caller.
185  friend NURBSPatch *Revolve3D(NURBSPatch &patch, double n[], double ang,
186  int times);
187 };
188 
189 
190 #ifdef MFEM_USE_MPI
191 class ParNURBSExtension;
192 #endif
193 
194 class NURBSPatchMap;
195 
196 
198 {
199 #ifdef MFEM_USE_MPI
200  friend class ParNURBSExtension;
201 #endif
202  friend class NURBSPatchMap;
203 
204 protected:
205  int mOrder; // see GetOrder() for description
208  // global entity counts
210  // local entity counts
213 
214  Array<int> activeVert; // activeVert[glob_vert] = loc_vert or -1
217  Array<int> activeDof; // activeDof[glob_dof] = loc_dof + 1 or 0
218 
220  int own_topo;
222  /** Set of knotvectors containing unique KnotVectors only */
224  /** Comprehensive set of knotvectors. This set contains a KnotVector for
225  every edge.*/
228 
229  // periodic BC info:
230  // - dof 2 dof map
231  // - master and slave boundary indices
235 
236  // global offsets, meshOffsets == meshVertexOffsets
241 
242  // global offsets, spaceOffsets == dofOffsets
247 
249 
252  Array2D<int> el_to_IJK; // IJK are "knot-span" indices!
253  Array2D<int> bel_to_IJK; // they are NOT element indices!
254 
255  std::vector<Array<int>> patch_to_el;
256  std::vector<Array<int>> patch_to_bel;
257 
259 
260  inline int KnotInd(int edge) const;
261  /// @note The returned object should NOT be deleted by the caller.
262  inline KnotVector *KnotVec(int edge);
263  inline const KnotVector *KnotVec(int edge) const;
264  inline const KnotVector *KnotVec(int edge, int oedge, int *okv) const;
265 
266  void CheckPatches();
267  void CheckBdrPatches();
268 
269  /** Checks the direction of the knotvectors in the patch based on
270  the patch orientation for patch @a p returns the direction of
271  the Knotvectors in @a kvdir.*/
272  void CheckKVDirection(int p, Array <int> &kvdir);
273  /** Creates the comprehensive set of KnotVectors. They are the same for 1D. */
274  void CreateComprehensiveKV();
275  /** Updates the unique set of KnotVectors */
276  void UpdateUniqueKV();
277 
278  /** Checks if the comprehensive array of KnotVectors agrees with
279  the reduced set of KnotVectors. Returns false if it finds
280  a difference. */
281  bool ConsistentKVSets();
282 
283  void GetPatchKnotVectors (int p, Array<KnotVector *> &kv);
285 
286  void SetOrderFromOrders();
288 
289  // periodic BC helper functions
290  void InitDofMap();
291  void ConnectBoundaries();
292  void ConnectBoundaries1D(int bnd0, int bnd1);
293  void ConnectBoundaries2D(int bnd0, int bnd1);
294  void ConnectBoundaries3D(int bnd0, int bnd1);
295  int DofMap(int dof) const
296  {
297  return (d_to_d.Size() > 0 )? d_to_d[dof] : dof;
298  };
299 
300  // also count the global NumOfVertices and the global NumOfDofs
301  void GenerateOffsets();
302  // count the global NumOfElements
303  void CountElements();
304  // count the global NumOfBdrElements
305  void CountBdrElements();
306 
307  // generate the mesh elements
308  void Get1DElementTopo(Array<Element *> &elements) const;
309  void Get2DElementTopo(Array<Element *> &elements) const;
310  void Get3DElementTopo(Array<Element *> &elements) const;
311 
312  // generate the boundary mesh elements
313  void Get1DBdrElementTopo(Array<Element *> &boundary) const;
314  void Get2DBdrElementTopo(Array<Element *> &boundary) const;
315  void Get3DBdrElementTopo(Array<Element *> &boundary) const;
316 
317  // FE space generation functions
318 
319  // based on activeElem, count NumOfActiveDofs, generate el_dof,
320  // el_to_patch, el_to_IJK, activeDof map (global-to-local)
322 
323  // generate elem_to_global-dof table for the active elements
324  // define el_to_patch, el_to_IJK, activeDof (as bool)
328 
329  // call after GenerateElementDofTable
331 
332  // generate the bdr-elem_to_global-dof table for the active bdr. elements
333  // define bel_to_patch, bel_to_IJK
337 
338  // FE --> Patch translation functions
339  void GetPatchNets (const Vector &Nodes, int vdim);
340  void Get1DPatchNets(const Vector &Nodes, int vdim);
341  void Get2DPatchNets(const Vector &Nodes, int vdim);
342  void Get3DPatchNets(const Vector &Nodes, int vdim);
343 
344  // Patch --> FE translation functions
345  // Side effects: delete the patches, update the weights from the patches
346  void SetSolutionVector (Vector &Nodes, int vdim);
347  void Set1DSolutionVector(Vector &Nodes, int vdim);
348  void Set2DSolutionVector(Vector &Nodes, int vdim);
349  void Set3DSolutionVector(Vector &Nodes, int vdim);
350 
351  // determine activeVert, NumOfActiveVertices from the activeElem array
352  void GenerateActiveVertices();
353 
354  // determine activeBdrElem, NumOfActiveBdrElems
355  void GenerateActiveBdrElems();
356 
357  void MergeWeights(Mesh *mesh_array[], int num_pieces);
358 
359  void SetPatchToElements();
360  void SetPatchToBdrElements();
361 
362  // to be used by ParNURBSExtension constructor(s)
364 
365 public:
366  /// Copy constructor: deep copy
367  NURBSExtension(const NURBSExtension &orig);
368  /// Read-in a NURBSExtension
369  NURBSExtension(std::istream &input);
370  /** @brief Create a NURBSExtension with elevated order by repeating the
371  endpoints of the knot vectors and using uniform weights of 1. */
372  /** If a knot vector in @a parent already has order greater than or equal to
373  @a newOrder, it will be used unmodified. */
374  NURBSExtension(NURBSExtension *parent, int newOrder);
375  /** @brief Create a NURBSExtension with elevated knot vector orders (by
376  repeating the endpoints of the knot vectors and using uniform weights of
377  1) as given by the array @a newOrders. */
378  /** If a knot vector in @a parent already has order greater than or equal to
379  the corresponding entry in @a newOrder, it will be used unmodified. */
380  NURBSExtension(NURBSExtension *parent, const Array<int> &newOrders);
381  /// Construct a NURBSExtension by merging a partitioned NURBS mesh
382  NURBSExtension(Mesh *mesh_array[], int num_pieces);
383 
384  /// Copy assignment not supported
385  NURBSExtension& operator=(const NURBSExtension&) = delete;
386 
387  // Generate connections between boundaries, such as periodic BCs
389  const Array<int> &GetMaster() const { return master; };
390  Array<int> &GetMaster() { return master; };
391  const Array<int> &GetSlave() const { return slave; };
392  Array<int> &GetSlave() { return slave; };
393  void MergeGridFunctions(GridFunction *gf_array[], int num_pieces,
394  GridFunction &merged);
395 
396  /// Destroy a NURBSExtension
397  virtual ~NURBSExtension();
398 
399  // Print functions
400  void Print(std::ostream &out) const;
401  void PrintCharacteristics(std::ostream &out) const;
402  void PrintFunctions(const char *filename, int samples=11) const;
403 
404  // Meta data functions
405  int Dimension() const { return patchTopo->Dimension(); }
406  int GetNP() const { return patchTopo->GetNE(); }
407  int GetNBP() const { return patchTopo->GetNBE(); }
408 
409  /// Read-only access to the orders of all knot vectors.
410  const Array<int> &GetOrders() const { return mOrders; }
411  /** @brief If all orders are identical, return that number. Otherwise, return
412  NURBSFECollection::VariableOrder. */
413  int GetOrder() const { return mOrder; }
414 
415  int GetNKV() const { return NumOfKnotVectors; }
416 
417  int GetGNV() const { return NumOfVertices; }
418  int GetNV() const { return NumOfActiveVertices; }
419  int GetGNE() const { return NumOfElements; }
420  int GetNE() const { return NumOfActiveElems; }
421  int GetGNBE() const { return NumOfBdrElements; }
422  int GetNBE() const { return NumOfActiveBdrElems; }
423 
424  int GetNTotalDof() const { return NumOfDofs; }
425  int GetNDof() const { return NumOfActiveDofs; }
426 
427  /// Returns knotvectors in each dimension for patch @a p.
428  void GetPatchKnotVectors(int p, Array<const KnotVector *> &kv) const;
429 
431 
432  // Knotvector read-only access function
433  const KnotVector *GetKnotVector(int i) const { return knotVectors[i]; }
434 
435  // Mesh generation functions
436  void GetElementTopo (Array<Element *> &elements) const;
437  void GetBdrElementTopo(Array<Element *> &boundary) const;
438 
439  bool HavePatches() const { return (patches.Size() != 0); }
440 
441  /// @note The returned object should NOT be deleted by the caller.
443  /// @note The returned object should NOT be deleted by the caller.
445 
446  void GetVertexLocalToGlobal(Array<int> &lvert_vert);
447  void GetElementLocalToGlobal(Array<int> &lelem_elem);
448 
449  // Set the attribute for patch @a i, which is set to all elements in the
450  // patch.
451  void SetPatchAttribute(int i, int attr) { patchTopo->SetAttribute(i, attr); }
452 
453  // Get the attribute for patch @a i, which is set to all elements in the
454  // patch.
455  int GetPatchAttribute(int i) const { return patchTopo->GetAttribute(i); }
456 
457  // Set the attribute for patch boundary element @a i, which is set to all
458  // boundary elements in the patch.
459  void SetPatchBdrAttribute(int i, int attr)
460  { patchTopo->SetBdrAttribute(i, attr); }
461  // Get the attribute for patch boundary element @a i, which is set to all
462  // boundary elements in the patch.
463  int GetPatchBdrAttribute(int i) const
464  { return patchTopo->GetBdrAttribute(i); }
465 
466  // Load functions
467  void LoadFE(int i, const FiniteElement *FE) const;
468  void LoadBE(int i, const FiniteElement *BE) const;
469 
470  const Vector &GetWeights() const { return weights; }
471  Vector &GetWeights() { return weights; }
472 
473  // Translation functions: from FE coordinates to IJK patch
474  // format and vice versa
475  void ConvertToPatches(const Vector &Nodes);
476  void SetKnotsFromPatches();
477  void SetCoordsFromPatches(Vector &Nodes);
478 
479  // Read a GridFunction written patch-by-patch, e.g. with PrintSolution().
480  void LoadSolution(std::istream &input, GridFunction &sol) const;
481  // Write a GridFunction patch-by-patch.
482  void PrintSolution(const GridFunction &sol, std::ostream &out) const;
483 
484  // Refinement methods
485  // new_degree = max(old_degree, min(old_degree + rel_degree, degree))
486  void DegreeElevate(int rel_degree, int degree = 16);
487  void UniformRefinement();
489  void KnotInsert(Array<Vector *> &kv);
490 
491  /// Returns the index of the patch containing element @a elem.
492  int GetElementPatch(int elem) const { return el_to_patch[elem]; }
493 
494  /** Returns the Cartesian indices (i,j) in 2D or (i,j,k) in 3D of element
495  @a elem, in the knot-span tensor product ordering for its patch. */
496  void GetElementIJK(int elem, Array<int> & ijk);
497 
498  // Returns the degrees of freedom on the patch, in Cartesian order.
499  void GetPatchDofs(const int patch, Array<int> &dofs);
500 
501  const Array<int>& GetPatchElements(int patch);
502  const Array<int>& GetPatchBdrElements(int patch);
503 };
504 
505 
506 #ifdef MFEM_USE_MPI
508 {
509 private:
510  int *partitioning;
511 
512  Table *GetGlobalElementDofTable();
513  Table *Get1DGlobalElementDofTable();
514  Table *Get2DGlobalElementDofTable();
515  Table *Get3DGlobalElementDofTable();
516 
517  void SetActive(const int *partitioning, const Array<bool> &active_bel);
518  void BuildGroups(const int *partitioning, const Table &elem_dof);
519 
520 public:
522 
524 
526 
527  ParNURBSExtension(MPI_Comm comm, NURBSExtension *parent, int *partitioning,
528  const Array<bool> &active_bel);
529 
530  // Create a parallel version of 'parent' with partitioning as in
531  // 'par_parent'; the 'parent' object is destroyed.
532  // The 'parent' can be either a local NURBSExtension or a global one.
534  const ParNURBSExtension *par_parent);
535 
536  virtual ~ParNURBSExtension() { delete [] partitioning; }
537 };
538 #endif
539 
540 
542 {
543 private:
544  const NURBSExtension *Ext;
545 
546  int I, J, K, pOffset, opatch;
547  Array<int> verts, edges, faces, oedge, oface;
548 
549  inline static int F(const int n, const int N)
550  { return (n < 0) ? 0 : ((n >= N) ? 2 : 1); }
551 
552  inline static int Or1D(const int n, const int N, const int Or)
553  { return (Or > 0) ? n : (N - 1 - n); }
554 
555  inline static int Or2D(const int n1, const int n2,
556  const int N1, const int N2, const int Or);
557 
558  // also set verts, edges, faces, orientations etc
559  void GetPatchKnotVectors (int p, const KnotVector *kv[]);
560  void GetBdrPatchKnotVectors(int p, const KnotVector *kv[], int *okv);
561 
562 public:
563  NURBSPatchMap(const NURBSExtension *ext) { Ext = ext; }
564 
565  int nx() { return I + 1; }
566  int ny() { return J + 1; }
567  int nz() { return K + 1; }
568 
569  void SetPatchVertexMap(int p, const KnotVector *kv[]);
570  void SetPatchDofMap (int p, const KnotVector *kv[]);
571 
572  void SetBdrPatchVertexMap(int p, const KnotVector *kv[], int *okv);
573  void SetBdrPatchDofMap (int p, const KnotVector *kv[], int *okv);
574 
575  inline int operator()(const int i) const;
576  inline int operator[](const int i) const { return (*this)(i); }
577 
578  inline int operator()(const int i, const int j) const;
579 
580  inline int operator()(const int i, const int j, const int k) const;
581 };
582 
583 
584 // Inline function implementations
585 
586 inline double &NURBSPatch::slice(int i, int j)
587 {
588 #ifdef MFEM_DEBUG
589  if (data == 0 || i < 0 || i >= nd || j < 0 || j > ls)
590  {
591  mfem_error("NURBSPatch::slice()");
592  }
593 #endif
594  return data[j%sd + sd*(i + (j/sd)*nd)];
595 }
596 
597 inline const double &NURBSPatch::slice(int i, int j) const
598 {
599 #ifdef MFEM_DEBUG
600  if (data == 0 || i < 0 || i >= nd || j < 0 || j > ls)
601  {
602  mfem_error("NURBSPatch::slice()");
603  }
604 #endif
605  return data[j%sd + sd*(i + (j/sd)*nd)];
606 }
607 
608 
609 inline double &NURBSPatch::operator()(int i, int l)
610 {
611 #ifdef MFEM_DEBUG
612  if (data == 0 || i < 0 || i >= ni || nj > 0 || nk > 0 ||
613  l < 0 || l >= Dim)
614  {
615  mfem_error("NURBSPatch::operator() 1D");
616  }
617 #endif
618 
619  return data[i*Dim+l];
620 }
621 
622 inline const double &NURBSPatch::operator()(int i, int l) const
623 {
624 #ifdef MFEM_DEBUG
625  if (data == 0 || i < 0 || i >= ni || nj > 0 || nk > 0 ||
626  l < 0 || l >= Dim)
627  {
628  mfem_error("NURBSPatch::operator() const 1D");
629  }
630 #endif
631 
632  return data[i*Dim+l];
633 }
634 
635 inline double &NURBSPatch::operator()(int i, int j, int l)
636 {
637 #ifdef MFEM_DEBUG
638  if (data == 0 || i < 0 || i >= ni || j < 0 || j >= nj || nk > 0 ||
639  l < 0 || l >= Dim)
640  {
641  mfem_error("NURBSPatch::operator() 2D");
642  }
643 #endif
644 
645  return data[(i+j*ni)*Dim+l];
646 }
647 
648 inline const double &NURBSPatch::operator()(int i, int j, int l) const
649 {
650 #ifdef MFEM_DEBUG
651  if (data == 0 || i < 0 || i >= ni || j < 0 || j >= nj || nk > 0 ||
652  l < 0 || l >= Dim)
653  {
654  mfem_error("NURBSPatch::operator() const 2D");
655  }
656 #endif
657 
658  return data[(i+j*ni)*Dim+l];
659 }
660 
661 inline double &NURBSPatch::operator()(int i, int j, int k, int l)
662 {
663 #ifdef MFEM_DEBUG
664  if (data == 0 || i < 0 || i >= ni || j < 0 || j >= nj || k < 0 ||
665  k >= nk || l < 0 || l >= Dim)
666  {
667  mfem_error("NURBSPatch::operator() 3D");
668  }
669 #endif
670 
671  return data[(i+(j+k*nj)*ni)*Dim+l];
672 }
673 
674 inline const double &NURBSPatch::operator()(int i, int j, int k, int l) const
675 {
676 #ifdef MFEM_DEBUG
677  if (data == 0 || i < 0 || i >= ni || j < 0 || j >= nj || k < 0 ||
678  k >= nk || l < 0 || l >= Dim)
679  {
680  mfem_error("NURBSPatch::operator() const 3D");
681  }
682 #endif
683 
684  return data[(i+(j+k*nj)*ni)*Dim+l];
685 }
686 
687 
688 inline int NURBSExtension::KnotInd(int edge) const
689 {
690  int kv = edge_to_knot[edge];
691  return (kv >= 0) ? kv : (-1-kv);
692 }
693 
695 {
696  return knotVectors[KnotInd(edge)];
697 }
698 
699 inline const KnotVector *NURBSExtension::KnotVec(int edge) const
700 {
701  return knotVectors[KnotInd(edge)];
702 }
703 
704 inline const KnotVector *NURBSExtension::KnotVec(int edge, int oedge, int *okv)
705 const
706 {
707  int kv = edge_to_knot[edge];
708  if (kv >= 0)
709  {
710  *okv = oedge;
711  return knotVectors[kv];
712  }
713  else
714  {
715  *okv = -oedge;
716  return knotVectors[-1-kv];
717  }
718 }
719 
720 
721 // static method
722 inline int NURBSPatchMap::Or2D(const int n1, const int n2,
723  const int N1, const int N2, const int Or)
724 {
725  // Needs testing
726  switch (Or)
727  {
728  case 0: return n1 + n2*N1;
729  case 1: return n2 + n1*N2;
730  case 2: return n2 + (N1 - 1 - n1)*N2;
731  case 3: return (N1 - 1 - n1) + n2*N1;
732  case 4: return (N1 - 1 - n1) + (N2 - 1 - n2)*N1;
733  case 5: return (N2 - 1 - n2) + (N1 - 1 - n1)*N2;
734  case 6: return (N2 - 1 - n2) + n1*N2;
735  case 7: return n1 + (N2 - 1 - n2)*N1;
736  }
737 #ifdef MFEM_DEBUG
738  mfem_error("NURBSPatchMap::Or2D");
739 #endif
740  return -1;
741 }
742 
743 inline int NURBSPatchMap::operator()(const int i) const
744 {
745  const int i1 = i - 1;
746  switch (F(i1, I))
747  {
748  case 0: return verts[0];
749  case 1: return pOffset + Or1D(i1, I, opatch);
750  case 2: return verts[1];
751  }
752 #ifdef MFEM_DEBUG
753  mfem_error("NURBSPatchMap::operator() const 1D");
754 #endif
755  return -1;
756 }
757 
758 inline int NURBSPatchMap::operator()(const int i, const int j) const
759 {
760  const int i1 = i - 1, j1 = j - 1;
761  switch (3*F(j1, J) + F(i1, I))
762  {
763  case 0: return verts[0];
764  case 1: return edges[0] + Or1D(i1, I, oedge[0]);
765  case 2: return verts[1];
766  case 3: return edges[3] + Or1D(j1, J, -oedge[3]);
767  case 4: return pOffset + Or2D(i1, j1, I, J, opatch);
768  case 5: return edges[1] + Or1D(j1, J, oedge[1]);
769  case 6: return verts[3];
770  case 7: return edges[2] + Or1D(i1, I, -oedge[2]);
771  case 8: return verts[2];
772  }
773 #ifdef MFEM_DEBUG
774  mfem_error("NURBSPatchMap::operator() const 2D");
775 #endif
776  return -1;
777 }
778 
779 inline int NURBSPatchMap::operator()(const int i, const int j, const int k)
780 const
781 {
782  // Needs testing
783  const int i1 = i - 1, j1 = j - 1, k1 = k - 1;
784  switch (3*(3*F(k1, K) + F(j1, J)) + F(i1, I))
785  {
786  case 0: return verts[0];
787  case 1: return edges[0] + Or1D(i1, I, oedge[0]);
788  case 2: return verts[1];
789  case 3: return edges[3] + Or1D(j1, J, oedge[3]);
790  case 4: return faces[0] + Or2D(i1, J - 1 - j1, I, J, oface[0]);
791  case 5: return edges[1] + Or1D(j1, J, oedge[1]);
792  case 6: return verts[3];
793  case 7: return edges[2] + Or1D(i1, I, oedge[2]);
794  case 8: return verts[2];
795  case 9: return edges[8] + Or1D(k1, K, oedge[8]);
796  case 10: return faces[1] + Or2D(i1, k1, I, K, oface[1]);
797  case 11: return edges[9] + Or1D(k1, K, oedge[9]);
798  case 12: return faces[4] + Or2D(J - 1 - j1, k1, J, K, oface[4]);
799  case 13: return pOffset + I*(J*k1 + j1) + i1;
800  case 14: return faces[2] + Or2D(j1, k1, J, K, oface[2]);
801  case 15: return edges[11] + Or1D(k1, K, oedge[11]);
802  case 16: return faces[3] + Or2D(I - 1 - i1, k1, I, K, oface[3]);
803  case 17: return edges[10] + Or1D(k1, K, oedge[10]);
804  case 18: return verts[4];
805  case 19: return edges[4] + Or1D(i1, I, oedge[4]);
806  case 20: return verts[5];
807  case 21: return edges[7] + Or1D(j1, J, oedge[7]);
808  case 22: return faces[5] + Or2D(i1, j1, I, J, oface[5]);
809  case 23: return edges[5] + Or1D(j1, J, oedge[5]);
810  case 24: return verts[7];
811  case 25: return edges[6] + Or1D(i1, I, oedge[6]);
812  case 26: return verts[6];
813  }
814 #ifdef MFEM_DEBUG
815  mfem_error("NURBSPatchMap::operator() const 3D");
816 #endif
817  return -1;
818 }
819 
820 }
821 
822 #endif
Array< KnotVector * > knotVectors
Definition: nurbs.hpp:223
Abstract class for all finite elements.
Definition: fe_base.hpp:233
Array< KnotVector * > kv
Definition: nurbs.hpp:107
int GetElementPatch(int elem) const
Returns the index of the patch containing element elem.
Definition: nurbs.hpp:492
Array< int > f_meshOffsets
Definition: nurbs.hpp:239
NURBSPatch(NURBSPatch *parent, int dir, int Order, int NCP)
Definition: nurbs.cpp:646
Array2D< int > bel_to_IJK
Definition: nurbs.hpp:253
void SetPatchAttribute(int i, int attr)
Definition: nurbs.hpp:451
void SetCoordsFromPatches(Vector &Nodes)
Definition: nurbs.cpp:3634
void Rotate2D(double angle)
Definition: nurbs.cpp:1284
Array< KnotVector * > knotVectorsCompr
Definition: nurbs.hpp:226
Array< int > activeVert
Definition: nurbs.hpp:214
void ConnectBoundaries2D(int bnd0, int bnd1)
Definition: nurbs.cpp:2080
KnotVector * KnotVec(int edge)
Definition: nurbs.hpp:694
void swap(NURBSPatch *np)
Definition: nurbs.cpp:661
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
void DegreeElevate(int dir, int t)
Definition: nurbs.cpp:966
int NumOfElements
Definition: nurbs.hpp:38
void Generate2DElementDofTable()
Definition: nurbs.cpp:3234
void MergeGridFunctions(GridFunction *gf_array[], int num_pieces, GridFunction &merged)
Definition: nurbs.cpp:2317
Array< int > ldof_group
Definition: nurbs.hpp:523
bool HavePatches() const
Definition: nurbs.hpp:439
int GetNV() const
Definition: nurbs.hpp:418
Array< int > master
Definition: nurbs.hpp:233
void PrintSolution(const GridFunction &sol, std::ostream &out) const
Definition: nurbs.cpp:3720
static void Get2DRotationMatrix(double angle, DenseMatrix &T)
Definition: nurbs.cpp:1272
virtual ~ParNURBSExtension()
Definition: nurbs.hpp:536
void LoadFE(int i, const FiniteElement *FE) const
Definition: nurbs.cpp:3579
void Generate3DElementDofTable()
Definition: nurbs.cpp:3288
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
NURBSExtension & operator=(const NURBSExtension &)=delete
Copy assignment not supported.
void KnotInsert(int dir, const KnotVector &knot)
Definition: nurbs.cpp:817
void MergeWeights(Mesh *mesh_array[], int num_pieces)
Definition: nurbs.cpp:2292
int GetNKV() const
Definition: nurbs.hpp:153
void Print(std::ostream &out) const
Definition: nurbs.cpp:126
NURBSPatchMap(const NURBSExtension *ext)
Definition: nurbs.hpp:563
void CalcDShape(Vector &grad, int i, double xi) const
Definition: nurbs.cpp:188
void UniformRefinement()
Definition: nurbs.cpp:3774
const Array< int > & GetSlave() const
Definition: nurbs.hpp:391
int DofMap(int dof) const
Definition: nurbs.hpp:295
void PrintCharacteristics(std::ostream &out) const
Definition: nurbs.cpp:1943
void Get2DBdrElementTopo(Array< Element *> &boundary) const
Definition: nurbs.cpp:3086
double & slice(int i, int j)
Definition: nurbs.hpp:586
void Rotate3D(double normal[], double angle)
Definition: nurbs.cpp:1349
Array< int > e_meshOffsets
Definition: nurbs.hpp:238
int Size() const
Returns the size of the vector.
Definition: vector.hpp:197
void GetElements()
Count the number of elements.
Definition: nurbs.cpp:101
void SetOrdersFromKnotVectors()
Definition: nurbs.cpp:2772
int GetOrder() const
Definition: nurbs.hpp:52
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
void GetPatchDofs(const int patch, Array< int > &dofs)
Definition: nurbs.cpp:3354
void Generate1DBdrElementDofTable()
Definition: nurbs.cpp:3428
int Dimension() const
Definition: nurbs.hpp:405
Array< int > & GetSlave()
Definition: nurbs.hpp:392
void Generate3DBdrElementDofTable()
Definition: nurbs.cpp:3501
int GetNE() const
Definition: nurbs.hpp:420
int NumOfControlPoints
Definition: nurbs.hpp:38
Array< int > e_spaceOffsets
Definition: nurbs.hpp:244
void GetBdrElementTopo(Array< Element *> &boundary) const
Definition: nurbs.cpp:3045
const Vector & GetWeights() const
Definition: nurbs.hpp:470
void GenerateBdrElementDofTable()
Definition: nurbs.cpp:3403
void CalcD2Shape(Vector &grad2, int i, double xi) const
Definition: nurbs.hpp:68
std::vector< Array< int > > patch_to_el
Definition: nurbs.hpp:255
void GetElementIJK(int elem, Array< int > &ijk)
Definition: nurbs.cpp:4096
int GetNBP() const
Definition: nurbs.hpp:407
void LoadSolution(std::istream &input, GridFunction &sol) const
Definition: nurbs.cpp:3683
void GetElementLocalToGlobal(Array< int > &lelem_elem)
Definition: nurbs.cpp:3569
void Get2DElementTopo(Array< Element *> &elements) const
Definition: nurbs.cpp:2964
void KnotInsert(Array< KnotVector *> &kv)
Definition: nurbs.cpp:3782
void Set1DSolutionVector(Vector &Nodes, int vdim)
Definition: nurbs.cpp:4010
void SetPatchToBdrElements()
Definition: nurbs.cpp:4113
void init(int dim_)
Definition: nurbs.cpp:504
void ConnectBoundaries1D(int bnd0, int bnd1)
Definition: nurbs.cpp:2066
Array< int > edge_to_knot
Definition: nurbs.hpp:221
void GenerateActiveVertices()
Definition: nurbs.cpp:2197
Array< bool > activeElem
Definition: nurbs.hpp:215
void Print(std::ostream &out) const
Definition: nurbs.cpp:1908
void GetVertexLocalToGlobal(Array< int > &lvert_vert)
Definition: nurbs.cpp:3559
void Set2DSolutionVector(Vector &Nodes, int vdim)
Definition: nurbs.cpp:4037
Vector knot
Definition: nurbs.hpp:37
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:1089
friend NURBSPatch * Revolve3D(NURBSPatch &patch, double n[], double ang, int times)
Definition: nurbs.cpp:1445
Array< int > el_to_patch
Definition: nurbs.hpp:250
Array< int > & GetMaster()
Definition: nurbs.hpp:390
int GetGNBE() const
Definition: nurbs.hpp:421
void SwapDirections(int dir1, int dir2)
Definition: nurbs.cpp:1230
KnotVector()
Create KnotVector.
Definition: nurbs.hpp:42
const Array< int > & GetMaster() const
Definition: nurbs.hpp:389
static const int MaxOrder
Definition: nurbs.hpp:35
Array< int > v_spaceOffsets
Definition: nurbs.hpp:243
void GetElementTopo(Array< Element *> &elements) const
Definition: nurbs.cpp:2916
void FindMaxima(Array< int > &ks, Vector &xi, Vector &u)
Definition: nurbs.cpp:346
int GetAttribute(int i) const
Return the attribute of element i.
Definition: mesh.hpp:1190
void Generate2DBdrElementDofTable()
Definition: nurbs.cpp:3458
int KnotInd(int edge) const
Definition: nurbs.hpp:688
void Generate1DElementDofTable()
Definition: nurbs.cpp:3191
void GetBdrPatchKnotVectors(int p, Array< KnotVector *> &kv)
Definition: nurbs.cpp:2717
int GetGNE() const
Definition: nurbs.hpp:419
void ConnectBoundaries3D(int bnd0, int bnd1)
Definition: nurbs.cpp:2126
ParNURBSExtension(const ParNURBSExtension &orig)
Definition: nurbs.cpp:4139
Array< int > d_to_d
Definition: nurbs.hpp:232
void Set3DSolutionVector(Vector &Nodes, int vdim)
Definition: nurbs.cpp:4065
void GenerateActiveBdrElems()
Definition: nurbs.cpp:2270
double * data
Definition: nurbs.hpp:105
int GetOrder() const
If all orders are identical, return that number. Otherwise, return NURBSFECollection::VariableOrder.
Definition: nurbs.hpp:413
void Get3DPatchNets(const Vector &Nodes, int vdim)
Definition: nurbs.cpp:3964
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:154
std::vector< Array< int > > patch_to_bel
Definition: nurbs.hpp:256
const Array< int > & GetOrders() const
Read-only access to the orders of all knot vectors.
Definition: nurbs.hpp:410
void Get2DPatchNets(const Vector &Nodes, int vdim)
Definition: nurbs.cpp:3937
Vector & GetWeights()
Definition: nurbs.hpp:471
int GetPatchAttribute(int i) const
Definition: nurbs.hpp:455
virtual ~NURBSExtension()
Destroy a NURBSExtension.
Definition: nurbs.cpp:1879
void GenerateOffsets()
Definition: nurbs.cpp:2782
int GetNCP() const
Definition: nurbs.hpp:51
void UniformRefinement(Vector &newknots) const
Definition: nurbs.cpp:87
Array< int > activeDof
Definition: nurbs.hpp:217
friend NURBSPatch * Interpolate(NURBSPatch &p1, NURBSPatch &p2)
Definition: nurbs.cpp:1398
void CalcDnShape(Vector &gradn, int n, int i, double xi) const
Definition: nurbs.cpp:245
static void Get3DRotationMatrix(double n[], double angle, double r, DenseMatrix &T)
Definition: nurbs.cpp:1310
void SetKnotsFromPatches()
Definition: nurbs.cpp:3642
int GetNTotalDof() const
Definition: nurbs.hpp:424
void UniformRefinement()
Definition: nurbs.cpp:799
void FindInterpolant(Array< Vector *> &x)
Definition: nurbs.cpp:409
Array< int > f_spaceOffsets
Definition: nurbs.hpp:245
Array< bool > activeBdrElem
Definition: nurbs.hpp:216
int operator()(const int i) const
Definition: nurbs.hpp:743
KnotVector(const KnotVector &kv)
Definition: nurbs.hpp:45
const Array< int > & GetPatchBdrElements(int patch)
Definition: nurbs.cpp:4131
void Get3DElementTopo(Array< Element *> &elements) const
Definition: nurbs.cpp:3000
void Get3DBdrElementTopo(Array< Element *> &boundary) const
Definition: nurbs.cpp:3117
int findKnotSpan(double u) const
Definition: nurbs.cpp:442
int GetPatchBdrAttribute(int i) const
Definition: nurbs.hpp:463
void SetPatchToElements()
Definition: nurbs.cpp:4102
KnotVector * GetKV(int i)
Definition: nurbs.hpp:155
Array< int > slave
Definition: nurbs.hpp:234
int GetBdrAttribute(int i) const
Return the attribute of boundary element i.
Definition: mesh.hpp:1196
int Size() const
Definition: nurbs.hpp:53
void Difference(const KnotVector &kv, Vector &diff) const
Definition: nurbs.cpp:471
int GetGNV() const
Definition: nurbs.hpp:417
int GetNE() const
Definition: nurbs.hpp:49
bool isElement(int i) const
Definition: nurbs.hpp:58
void LoadBE(int i, const FiniteElement *BE) const
Definition: nurbs.cpp:3600
int GetNC() const
Definition: nurbs.hpp:152
Array< int > v_meshOffsets
Definition: nurbs.hpp:237
int GetNKS() const
Definition: nurbs.hpp:50
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:66
void SetBdrPatchVertexMap(int p, const KnotVector *kv[], int *okv)
Definition: nurbs.cpp:4652
Array< int > p_meshOffsets
Definition: nurbs.hpp:240
void GetPatchKnotVectors(int p, Array< KnotVector *> &kv)
Definition: nurbs.cpp:2670
KnotVector & operator=(const KnotVector &kv)
Definition: nurbs.cpp:46
void SetPatchVertexMap(int p, const KnotVector *kv[])
Definition: nurbs.cpp:4589
void SetBdrAttribute(int i, int attr)
Set the attribute of boundary element i.
Definition: mesh.hpp:1199
const Array< int > & GetPatchElements(int patch)
Definition: nurbs.cpp:4124
Table * GetElementDofTable()
Definition: nurbs.hpp:442
Table * GetBdrElementDofTable()
Definition: nurbs.hpp:444
Array2D< int > el_to_IJK
Definition: nurbs.hpp:252
void SetOrderFromOrders()
Definition: nurbs.cpp:2758
void Get1DBdrElementTopo(Array< Element *> &boundary) const
Definition: nurbs.cpp:3063
const double & operator[](int i) const
Definition: nurbs.hpp:97
void CountBdrElements()
Definition: nurbs.cpp:2896
void CreateComprehensiveKV()
Definition: nurbs.cpp:2475
void Rotate(double angle, double normal[]=NULL)
Rotate the NURBSPatch.
Definition: nurbs.cpp:1255
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1086
void CheckKVDirection(int p, Array< int > &kvdir)
Definition: nurbs.cpp:2407
void SetPatchBdrAttribute(int i, int attr)
Definition: nurbs.hpp:459
void GetPatchNets(const Vector &Nodes, int vdim)
Definition: nurbs.cpp:3896
const KnotVector * GetKnotVector(int i) const
Definition: nurbs.hpp:433
int MakeUniformDegree(int degree=-1)
Definition: nurbs.cpp:1375
Array< int > p_spaceOffsets
Definition: nurbs.hpp:246
void SetSolutionVector(Vector &Nodes, int vdim)
Definition: nurbs.cpp:3994
void SetBdrPatchDofMap(int p, const KnotVector *kv[], int *okv)
Definition: nurbs.cpp:4685
void ConnectBoundaries()
Definition: nurbs.cpp:2000
double & operator()(int i, int j)
Definition: nurbs.hpp:609
void CalcShape(Vector &shape, int i, double xi) const
Definition: nurbs.cpp:161
void GenerateElementDofTable()
Definition: nurbs.cpp:3155
Array< NURBSPatch * > patches
Definition: nurbs.hpp:258
~KnotVector()
Destroys KnotVector.
Definition: nurbs.hpp:94
void Get1DElementTopo(Array< Element *> &elements) const
Definition: nurbs.cpp:2934
NURBSPatch & operator=(const NURBSPatch &)=delete
int GetNKV() const
Definition: nurbs.hpp:415
int operator[](const int i) const
Definition: nurbs.hpp:576
Array< int > mOrders
Definition: nurbs.hpp:206
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
RefCoord t[3]
void SetAttribute(int i, int attr)
Set the attribute of element i.
Definition: mesh.hpp:1193
Vector data type.
Definition: vector.hpp:58
void PrintFunctions(const char *filename, int samples=11) const
Definition: nurbs.cpp:1973
double u(const Vector &xvec)
Definition: lor_mms.hpp:22
bool ConsistentKVSets()
Definition: nurbs.cpp:2603
void Get1DPatchNets(const Vector &Nodes, int vdim)
Definition: nurbs.cpp:3912
void DegreeElevate(int rel_degree, int degree=16)
Definition: nurbs.cpp:3758
void CheckBdrPatches()
Definition: nurbs.cpp:2379
void PrintFunctions(std::ostream &out, int samples=11) const
Definition: nurbs.cpp:133
int GetNBE() const
Definition: nurbs.hpp:422
double getKnotLocation(double xi, int ni) const
Definition: nurbs.hpp:60
double & operator[](int i)
Definition: nurbs.hpp:96
void FlipDirection(int dir)
Definition: nurbs.cpp:1218
void Print(std::ostream &out) const
Definition: nurbs.cpp:700
GroupTopology gtopo
Definition: nurbs.hpp:521
Array< int > bel_to_patch
Definition: nurbs.hpp:251
void SetPatchDofMap(int p, const KnotVector *kv[])
Definition: nurbs.cpp:4621
int GetNP() const
Definition: nurbs.hpp:406
void ConvertToPatches(const Vector &Nodes)
Definition: nurbs.cpp:3623
int GetNDof() const
Definition: nurbs.hpp:425
int SetLoopDirection(int dir)
Definition: nurbs.cpp:724
KnotVector * DegreeElevate(int t) const
Definition: nurbs.cpp:58