MFEM  v4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
nurbs.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 #include "mesh_headers.hpp"
13 #include "../fem/fem.hpp"
14 #include "../general/text.hpp"
15 
16 #include <algorithm>
17 #if defined(_MSC_VER) && (_MSC_VER < 1800)
18 #include <float.h>
19 #define copysign _copysign
20 #endif
21 
22 namespace mfem
23 {
24 
25 using namespace std;
26 
27 const int KnotVector::MaxOrder = 10;
28 
29 KnotVector::KnotVector(std::istream &input)
30 {
31  input >> Order >> NumOfControlPoints;
32  knot.Load(input, NumOfControlPoints + Order + 1);
33  GetElements();
34 }
35 
36 KnotVector::KnotVector(int Order_, int NCP)
37 {
38  Order = Order_;
39  NumOfControlPoints = NCP;
40  knot.SetSize(NumOfControlPoints + Order + 1);
41 
42  knot = -1.;
43 }
44 
46 {
47  Order = kv.Order;
48  NumOfControlPoints = kv.NumOfControlPoints;
49  NumOfElements = kv.NumOfElements;
50  knot = kv.knot;
51  // alternatively, re-compute NumOfElements
52  // GetElements();
53 
54  return *this;
55 }
56 
58 {
59  if (t < 0)
60  {
61  mfem_error("KnotVector::DegreeElevate :\n"
62  " Parent KnotVector order higher than child");
63  }
64 
65  const int nOrder = Order + t;
66  KnotVector *newkv = new KnotVector(nOrder, GetNCP() + t);
67 
68  for (int i = 0; i <= nOrder; i++)
69  {
70  (*newkv)[i] = knot(0);
71  }
72  for (int i = nOrder + 1; i < newkv->GetNCP(); i++)
73  {
74  (*newkv)[i] = knot(i - t);
75  }
76  for (int i = 0; i <= nOrder; i++)
77  {
78  (*newkv)[newkv->GetNCP() + i] = knot(knot.Size()-1);
79  }
80 
81  newkv->GetElements();
82 
83  return newkv;
84 }
85 
87 {
88  newknots.SetSize(NumOfElements);
89  int j = 0;
90  for (int i = 0; i < knot.Size()-1; i++)
91  {
92  if (knot(i) != knot(i+1))
93  {
94  newknots(j) = 0.5*(knot(i) + knot(i+1));
95  j++;
96  }
97  }
98 }
99 
101 {
102  NumOfElements = 0;
103  for (int i = Order; i < NumOfControlPoints; i++)
104  {
105  if (knot(i) != knot(i+1))
106  {
107  NumOfElements++;
108  }
109  }
110 }
111 
113 {
114  double apb = knot(0) + knot(knot.Size()-1);
115 
116  int ns = (NumOfControlPoints - Order)/2;
117  for (int i = 1; i <= ns; i++)
118  {
119  double tmp = apb - knot(Order + i);
120  knot(Order + i) = apb - knot(NumOfControlPoints - i);
121  knot(NumOfControlPoints - i) = tmp;
122  }
123 }
124 
125 void KnotVector::Print(std::ostream &out) const
126 {
127  out << Order << ' ' << NumOfControlPoints << ' ';
128  knot.Print(out, knot.Size());
129 }
130 
131 // Routine from "The NURBS book" - 2nd ed - Piegl and Tiller
132 void KnotVector::CalcShape(Vector &shape, int i, double xi) const
133 {
134  MFEM_ASSERT(Order <= MaxOrder, "Order > MaxOrder!");
135 
136  int p = Order;
137  int ip = (i >= 0) ? (i + p) : (-1 - i + p);
138  double u = getKnotLocation((i >= 0) ? xi : 1. - xi, ip), saved, tmp;
139  double left[MaxOrder+1], right[MaxOrder+1];
140 
141  shape(0) = 1.;
142  for (int j = 1; j <= p; ++j)
143  {
144  left[j] = u - knot(ip+1-j);
145  right[j] = knot(ip+j) - u;
146  saved = 0.;
147  for (int r = 0; r < j; ++r)
148  {
149  tmp = shape(r)/(right[r+1] + left[j-r]);
150  shape(r) = saved + right[r+1]*tmp;
151  saved = left[j-r]*tmp;
152  }
153  shape(j) = saved;
154  }
155 }
156 
157 // Routine from "The NURBS book" - 2nd ed - Piegl and Tiller
158 void KnotVector::CalcDShape(Vector &grad, int i, double xi) const
159 {
160  int p = Order, rk, pk;
161  int ip = (i >= 0) ? (i + p) : (-1 - i + p);
162  double u = getKnotLocation((i >= 0) ? xi : 1. - xi, ip), temp, saved, d;
163  double ndu[MaxOrder+1][MaxOrder+1], left[MaxOrder+1], right[MaxOrder+1];
164 
165 #ifdef MFEM_DEBUG
166  if (p > MaxOrder)
167  {
168  mfem_error("KnotVector::CalcDShape : Order > MaxOrder!");
169  }
170 #endif
171 
172  ndu[0][0] = 1.0;
173  for (int j = 1; j <= p; j++)
174  {
175  left[j] = u - knot(ip-j+1);
176  right[j] = knot(ip+j) - u;
177  saved = 0.0;
178  for (int r = 0; r < j; r++)
179  {
180  ndu[j][r] = right[r+1] + left[j-r];
181  temp = ndu[r][j-1]/ndu[j][r];
182  ndu[r][j] = saved + right[r+1]*temp;
183  saved = left[j-r]*temp;
184  }
185  ndu[j][j] = saved;
186  }
187 
188  for (int r = 0; r <= p; ++r)
189  {
190  d = 0.0;
191  rk = r-1;
192  pk = p-1;
193  if (r >= 1)
194  {
195  d = ndu[rk][pk]/ndu[p][rk];
196  }
197  if (r <= pk)
198  {
199  d -= ndu[r][pk]/ndu[p][r];
200  }
201  grad(r) = d;
202  }
203 
204  if (i >= 0)
205  {
206  grad *= p*(knot(ip+1) - knot(ip));
207  }
208  else
209  {
210  grad *= p*(knot(ip) - knot(ip+1));
211  }
212 }
213 
214 int KnotVector::findKnotSpan(double u) const
215 {
216  int low, mid, high;
217 
218  if (u == knot(NumOfControlPoints+Order))
219  {
220  mid = NumOfControlPoints;
221  }
222  else
223  {
224  low = Order;
225  high = NumOfControlPoints + 1;
226  mid = (low + high)/2;
227  while ( (u < knot(mid-1)) || (u > knot(mid)) )
228  {
229  if (u < knot(mid-1))
230  {
231  high = mid;
232  }
233  else
234  {
235  low = mid;
236  }
237  mid = (low + high)/2;
238  }
239  }
240  return mid;
241 }
242 
243 void KnotVector::Difference(const KnotVector &kv, Vector &diff) const
244 {
245  if (Order != kv.GetOrder())
246  {
247  mfem_error("KnotVector::Difference :\n"
248  " Can not compare knot vectors with different orders!");
249  }
250 
251  int s = kv.Size() - Size();
252  if (s < 0)
253  {
254  kv.Difference(*this, diff);
255  return;
256  }
257 
258  diff.SetSize(s);
259 
260  s = 0;
261  int i = 0;
262  for (int j = 0; j < kv.Size(); j++)
263  {
264  if (knot(i) == kv[j])
265  {
266  i++;
267  }
268  else
269  {
270  diff(s) = kv[j];
271  s++;
272  }
273  }
274 }
275 
276 
277 void NURBSPatch::init(int dim_)
278 {
279  Dim = dim_;
280  sd = nd = -1;
281 
282  if (kv.Size() == 2)
283  {
284  ni = kv[0]->GetNCP();
285  nj = kv[1]->GetNCP();
286  nk = -1;
287 
288  data = new double[ni*nj*Dim];
289 
290 #ifdef MFEM_DEBUG
291  for (int i = 0; i < ni*nj*Dim; i++)
292  {
293  data[i] = -999.99;
294  }
295 #endif
296  }
297  else if (kv.Size() == 3)
298  {
299  ni = kv[0]->GetNCP();
300  nj = kv[1]->GetNCP();
301  nk = kv[2]->GetNCP();
302 
303  data = new double[ni*nj*nk*Dim];
304 
305 #ifdef MFEM_DEBUG
306  for (int i = 0; i < ni*nj*nk*Dim; i++)
307  {
308  data[i] = -999.99;
309  }
310 #endif
311  }
312  else
313  {
314  mfem_error("NURBSPatch::init : Wrond dimension of knotvectors!");
315  }
316 }
317 
318 NURBSPatch::NURBSPatch(const NURBSPatch &orig)
319  : ni(orig.ni), nj(orig.nj), nk(orig.nk), Dim(orig.Dim),
320  data(NULL), kv(orig.kv.Size()), sd(orig.sd), nd(orig.nd)
321 {
322  // Allocate and copy data:
323  const int data_size = Dim*ni*nj*((kv.Size() == 2) ? 1 : nk);
324  data = new double[data_size];
325  std::memcpy(data, orig.data, data_size*sizeof(double));
326 
327  // Copy the knot vectors:
328  for (int i = 0; i < kv.Size(); i++)
329  {
330  kv[i] = new KnotVector(*orig.kv[i]);
331  }
332 }
333 
334 NURBSPatch::NURBSPatch(std::istream &input)
335 {
336  int pdim, dim, size = 1;
337  string ident;
338 
339  input >> ws >> ident >> pdim; // knotvectors
340  kv.SetSize(pdim);
341  for (int i = 0; i < pdim; i++)
342  {
343  kv[i] = new KnotVector(input);
344  size *= kv[i]->GetNCP();
345  }
346 
347  input >> ws >> ident >> dim; // dimension
348  init(dim + 1);
349 
350  input >> ws >> ident; // controlpoints (homogeneous coordinates)
351  if (ident == "controlpoints" || ident == "controlpoints_homogeneous")
352  {
353  for (int j = 0, i = 0; i < size; i++)
354  for (int d = 0; d <= dim; d++, j++)
355  {
356  input >> data[j];
357  }
358  }
359  else // "controlpoints_cartesian" (Cartesian coordinates with weight)
360  {
361  for (int j = 0, i = 0; i < size; i++)
362  {
363  for (int d = 0; d <= dim; d++)
364  {
365  input >> data[j+d];
366  }
367  for (int d = 0; d < dim; d++)
368  {
369  data[j+d] *= data[j+dim];
370  }
371  j += (dim+1);
372  }
373  }
374 }
375 
376 NURBSPatch::NURBSPatch(const KnotVector *kv0, const KnotVector *kv1, int dim_)
377 {
378  kv.SetSize(2);
379  kv[0] = new KnotVector(*kv0);
380  kv[1] = new KnotVector(*kv1);
381  init(dim_);
382 }
383 
384 NURBSPatch::NURBSPatch(const KnotVector *kv0, const KnotVector *kv1,
385  const KnotVector *kv2, int dim_)
386 {
387  kv.SetSize(3);
388  kv[0] = new KnotVector(*kv0);
389  kv[1] = new KnotVector(*kv1);
390  kv[2] = new KnotVector(*kv2);
391  init(dim_);
392 }
393 
394 NURBSPatch::NURBSPatch(Array<const KnotVector *> &kv_, int dim_)
395 {
396  kv.SetSize(kv_.Size());
397  for (int i = 0; i < kv.Size(); i++)
398  {
399  kv[i] = new KnotVector(*kv_[i]);
400  }
401  init(dim_);
402 }
403 
404 NURBSPatch::NURBSPatch(NURBSPatch *parent, int dir, int Order, int NCP)
405 {
406  kv.SetSize(parent->kv.Size());
407  for (int i = 0; i < kv.Size(); i++)
408  if (i != dir)
409  {
410  kv[i] = new KnotVector(*parent->kv[i]);
411  }
412  else
413  {
414  kv[i] = new KnotVector(Order, NCP);
415  }
416  init(parent->Dim);
417 }
418 
419 void NURBSPatch::swap(NURBSPatch *np)
420 {
421  if (data != NULL)
422  {
423  delete [] data;
424  }
425 
426  for (int i = 0; i < kv.Size(); i++)
427  {
428  if (kv[i]) { delete kv[i]; }
429  }
430 
431  data = np->data;
432  np->kv.Copy(kv);
433 
434  ni = np->ni;
435  nj = np->nj;
436  nk = np->nk;
437  Dim = np->Dim;
438 
439  np->data = NULL;
440  np->kv.SetSize(0);
441 
442  delete np;
443 }
444 
445 NURBSPatch::~NURBSPatch()
446 {
447  if (data != NULL)
448  {
449  delete [] data;
450  }
451 
452  for (int i = 0; i < kv.Size(); i++)
453  {
454  if (kv[i]) { delete kv[i]; }
455  }
456 }
457 
458 void NURBSPatch::Print(std::ostream &out) const
459 {
460  int size = 1;
461 
462  out << "knotvectors\n" << kv.Size() << '\n';
463  for (int i = 0; i < kv.Size(); i++)
464  {
465  kv[i]->Print(out);
466  size *= kv[i]->GetNCP();
467  }
468 
469  out << "\ndimension\n" << Dim - 1
470  << "\n\ncontrolpoints\n";
471  for (int j = 0, i = 0; i < size; i++)
472  {
473  out << data[j++];
474  for (int d = 1; d < Dim; d++)
475  {
476  out << ' ' << data[j++];
477  }
478  out << '\n';
479  }
480 }
481 
482 int NURBSPatch::SetLoopDirection(int dir)
483 {
484  if (nk == -1)
485  {
486  if (dir == 0)
487  {
488  sd = Dim;
489  nd = ni;
490 
491  return nj*Dim;
492  }
493  else if (dir == 1)
494  {
495  sd = ni*Dim;
496  nd = nj;
497 
498  return ni*Dim;
499  }
500  else
501  {
502  mfem::err << "NURBSPatch::SetLoopDirection :\n"
503  " Direction error in 2D patch, dir = " << dir << '\n';
504  mfem_error();
505  }
506  }
507  else
508  {
509  if (dir == 0)
510  {
511  sd = Dim;
512  nd = ni;
513 
514  return nj*nk*Dim;
515  }
516  else if (dir == 1)
517  {
518  sd = ni*Dim;
519  nd = nj;
520 
521  return ni*nk*Dim;
522  }
523  else if (dir == 2)
524  {
525  sd = ni*nj*Dim;
526  nd = nk;
527 
528  return ni*nj*Dim;
529  }
530  else
531  {
532  mfem::err << "NURBSPatch::SetLoopDirection :\n"
533  " Direction error in 3D patch, dir = " << dir << '\n';
534  mfem_error();
535  }
536  }
537 
538  return -1;
539 }
540 
541 void NURBSPatch::UniformRefinement()
542 {
543  Vector newknots;
544  for (int dir = 0; dir < kv.Size(); dir++)
545  {
546  kv[dir]->UniformRefinement(newknots);
547  KnotInsert(dir, newknots);
548  }
549 }
550 
551 void NURBSPatch::KnotInsert(Array<KnotVector *> &newkv)
552 {
553  for (int dir = 0; dir < kv.Size(); dir++)
554  {
555  KnotInsert(dir, *newkv[dir]);
556  }
557 }
558 
559 void NURBSPatch::KnotInsert(int dir, const KnotVector &newkv)
560 {
561  if (dir >= kv.Size() || dir < 0)
562  {
563  mfem_error("NURBSPatch::KnotInsert : Incorrect direction!");
564  }
565 
566  int t = newkv.GetOrder() - kv[dir]->GetOrder();
567 
568  if (t > 0)
569  {
570  DegreeElevate(dir, t);
571  }
572  else if (t < 0)
573  {
574  mfem_error("NURBSPatch::KnotInsert : Incorrect order!");
575  }
576 
577  Vector diff;
578  GetKV(dir)->Difference(newkv, diff);
579  if (diff.Size() > 0)
580  {
581  KnotInsert(dir, diff);
582  }
583 }
584 
585 // Routine from "The NURBS book" - 2nd ed - Piegl and Tiller
586 void NURBSPatch::KnotInsert(int dir, const Vector &knot)
587 {
588  if (dir >= kv.Size() || dir < 0)
589  {
590  mfem_error("NURBSPatch::KnotInsert : Incorrect direction!");
591  }
592 
593  NURBSPatch &oldp = *this;
594  KnotVector &oldkv = *kv[dir];
595 
596  NURBSPatch *newpatch = new NURBSPatch(this, dir, oldkv.GetOrder(),
597  oldkv.GetNCP() + knot.Size());
598  NURBSPatch &newp = *newpatch;
599  KnotVector &newkv = *newp.GetKV(dir);
600 
601  int size = oldp.SetLoopDirection(dir);
602  if (size != newp.SetLoopDirection(dir))
603  {
604  mfem_error("NURBSPatch::KnotInsert : Size mismatch!");
605  }
606 
607  int rr = knot.Size() - 1;
608  int a = oldkv.findKnotSpan(knot(0)) - 1;
609  int b = oldkv.findKnotSpan(knot(rr)) - 1;
610  int pl = oldkv.GetOrder();
611  int ml = oldkv.GetNCP();
612 
613  for (int j = 0; j <= a; j++)
614  {
615  newkv[j] = oldkv[j];
616  }
617  for (int j = b+pl; j <= ml+pl; j++)
618  {
619  newkv[j+rr+1] = oldkv[j];
620  }
621  for (int k = 0; k <= (a-pl); k++)
622  {
623  for (int ll = 0; ll < size; ll++)
624  {
625  newp(k,ll) = oldp(k,ll);
626  }
627  }
628  for (int k = (b-1); k < ml; k++)
629  {
630  for (int ll = 0; ll < size; ll++)
631  {
632  newp(k+rr+1,ll) = oldp(k,ll);
633  }
634  }
635 
636  int i = b+pl-1;
637  int k = b+pl+rr;
638 
639  for (int j = rr; j >= 0; j--)
640  {
641  while ( (knot(j) <= oldkv[i]) && (i > a) )
642  {
643  newkv[k] = oldkv[i];
644  for (int ll = 0; ll < size; ll++)
645  {
646  newp(k-pl-1,ll) = oldp(i-pl-1,ll);
647  }
648 
649  k--;
650  i--;
651  }
652 
653  for (int ll = 0; ll < size; ll++)
654  {
655  newp(k-pl-1,ll) = newp(k-pl,ll);
656  }
657 
658  for (int l = 1; l <= pl; l++)
659  {
660  int ind = k-pl+l;
661  double alfa = newkv[k+l] - knot(j);
662  if (fabs(alfa) == 0.0)
663  {
664  for (int ll = 0; ll < size; ll++)
665  {
666  newp(ind-1,ll) = newp(ind,ll);
667  }
668  }
669  else
670  {
671  alfa = alfa/(newkv[k+l] - oldkv[i-pl+l]);
672  for (int ll = 0; ll < size; ll++)
673  {
674  newp(ind-1,ll) = alfa*newp(ind-1,ll) + (1.0-alfa)*newp(ind,ll);
675  }
676  }
677  }
678 
679  newkv[k] = knot(j);
680  k--;
681  }
682 
683  newkv.GetElements();
684 
685  swap(newpatch);
686 }
687 
688 void NURBSPatch::DegreeElevate(int t)
689 {
690  for (int dir = 0; dir < kv.Size(); dir++)
691  {
692  DegreeElevate(dir, t);
693  }
694 }
695 
696 // Routine from "The NURBS book" - 2nd ed - Piegl and Tiller
697 void NURBSPatch::DegreeElevate(int dir, int t)
698 {
699  if (dir >= kv.Size() || dir < 0)
700  {
701  mfem_error("NURBSPatch::DegreeElevate : Incorrect direction!");
702  }
703 
704  int i, j, k, kj, mpi, mul, mh, kind, cind, first, last;
705  int r, a, b, oldr, save, s, tr, lbz, rbz, l;
706  double inv, ua, ub, numer, alf, den, bet, gam;
707 
708  NURBSPatch &oldp = *this;
709  KnotVector &oldkv = *kv[dir];
710 
711  NURBSPatch *newpatch = new NURBSPatch(this, dir, oldkv.GetOrder() + t,
712  oldkv.GetNCP() + oldkv.GetNE()*t);
713  NURBSPatch &newp = *newpatch;
714  KnotVector &newkv = *newp.GetKV(dir);
715 
716  int size = oldp.SetLoopDirection(dir);
717  if (size != newp.SetLoopDirection(dir))
718  {
719  mfem_error("NURBSPatch::DegreeElevate : Size mismatch!");
720  }
721 
722  int p = oldkv.GetOrder();
723  int n = oldkv.GetNCP()-1;
724 
725  DenseMatrix bezalfs (p+t+1, p+1);
726  DenseMatrix bpts (p+1, size);
727  DenseMatrix ebpts (p+t+1, size);
728  DenseMatrix nextbpts(p-1, size);
729  Vector alphas (p-1);
730 
731  int m = n + p + 1;
732  int ph = p + t;
733  int ph2 = ph/2;
734 
735  {
736  Array2D<int> binom(ph+1, ph+1);
737  for (i = 0; i <= ph; i++)
738  {
739  binom(i,0) = binom(i,i) = 1;
740  for (j = 1; j < i; j++)
741  {
742  binom(i,j) = binom(i-1,j) + binom(i-1,j-1);
743  }
744  }
745 
746  bezalfs(0,0) = 1.0;
747  bezalfs(ph,p) = 1.0;
748 
749  for (i = 1; i <= ph2; i++)
750  {
751  inv = 1.0/binom(ph,i);
752  mpi = min(p,i);
753  for (j = max(0,i-t); j <= mpi; j++)
754  {
755  bezalfs(i,j) = inv*binom(p,j)*binom(t,i-j);
756  }
757  }
758  }
759 
760  for (i = ph2+1; i < ph; i++)
761  {
762  mpi = min(p,i);
763  for (j = max(0,i-t); j <= mpi; j++)
764  {
765  bezalfs(i,j) = bezalfs(ph-i,p-j);
766  }
767  }
768 
769  mh = ph;
770  kind = ph + 1;
771  r = -1;
772  a = p;
773  b = p + 1;
774  cind = 1;
775  ua = oldkv[0];
776  for (l = 0; l < size; l++)
777  {
778  newp(0,l) = oldp(0,l);
779  }
780  for (i = 0; i <= ph; i++)
781  {
782  newkv[i] = ua;
783  }
784 
785  for (i = 0; i <= p; i++)
786  {
787  for (l = 0; l < size; l++)
788  {
789  bpts(i,l) = oldp(i,l);
790  }
791  }
792 
793  while (b < m)
794  {
795  i = b;
796  while (b < m && oldkv[b] == oldkv[b+1]) { b++; }
797 
798  mul = b-i+1;
799 
800  mh = mh + mul + t;
801  ub = oldkv[b];
802  oldr = r;
803  r = p-mul;
804  if (oldr > 0) { lbz = (oldr+2)/2; }
805  else { lbz = 1; }
806 
807  if (r > 0) { rbz = ph-(r+1)/2; }
808  else { rbz = ph; }
809 
810  if (r > 0)
811  {
812  numer = ub - ua;
813  for (k = p ; k > mul; k--)
814  {
815  alphas[k-mul-1] = numer/(oldkv[a+k]-ua);
816  }
817 
818  for (j = 1; j <= r; j++)
819  {
820  save = r-j;
821  s = mul+j;
822  for (k = p; k >= s; k--)
823  {
824  for (l = 0; l < size; l++)
825  bpts(k,l) = (alphas[k-s]*bpts(k,l) +
826  (1.0-alphas[k-s])*bpts(k-1,l));
827  }
828  for (l = 0; l < size; l++)
829  {
830  nextbpts(save,l) = bpts(p,l);
831  }
832  }
833  }
834 
835  for (i = lbz; i <= ph; i++)
836  {
837  for (l = 0; l < size; l++)
838  {
839  ebpts(i,l) = 0.0;
840  }
841  mpi = min(p,i);
842  for (j = max(0,i-t); j <= mpi; j++)
843  {
844  for (l = 0; l < size; l++)
845  {
846  ebpts(i,l) += bezalfs(i,j)*bpts(j,l);
847  }
848  }
849  }
850 
851  if (oldr > 1)
852  {
853  first = kind-2;
854  last = kind;
855  den = ub-ua;
856  bet = (ub-newkv[kind-1])/den;
857 
858  for (tr = 1; tr < oldr; tr++)
859  {
860  i = first;
861  j = last;
862  kj = j-kind+1;
863  while (j-i > tr)
864  {
865  if (i < cind)
866  {
867  alf = (ub-newkv[i])/(ua-newkv[i]);
868  for (l = 0; l < size; l++)
869  {
870  newp(i,l) = alf*newp(i,l)-(1.0-alf)*newp(i-1,l);
871  }
872  }
873  if (j >= lbz)
874  {
875  if ((j-tr) <= (kind-ph+oldr))
876  {
877  gam = (ub-newkv[j-tr])/den;
878  for (l = 0; l < size; l++)
879  {
880  ebpts(kj,l) = gam*ebpts(kj,l) + (1.0-gam)*ebpts(kj+1,l);
881  }
882  }
883  else
884  {
885  for (l = 0; l < size; l++)
886  {
887  ebpts(kj,l) = bet*ebpts(kj,l) + (1.0-bet)*ebpts(kj+1,l);
888  }
889  }
890  }
891  i = i+1;
892  j = j-1;
893  kj = kj-1;
894  }
895  first--;
896  last++;
897  }
898  }
899 
900  if (a != p)
901  {
902  for (i = 0; i < (ph-oldr); i++)
903  {
904  newkv[kind] = ua;
905  kind = kind+1;
906  }
907  }
908  for (j = lbz; j <= rbz; j++)
909  {
910  for (l = 0; l < size; l++)
911  {
912  newp(cind,l) = ebpts(j,l);
913  }
914  cind = cind +1;
915  }
916 
917  if (b < m)
918  {
919  for (j = 0; j <r; j++)
920  for (l = 0; l < size; l++)
921  {
922  bpts(j,l) = nextbpts(j,l);
923  }
924 
925  for (j = r; j <= p; j++)
926  for (l = 0; l < size; l++)
927  {
928  bpts(j,l) = oldp(b-p+j,l);
929  }
930 
931  a = b;
932  b = b+1;
933  ua = ub;
934  }
935  else
936  {
937  for (i = 0; i <= ph; i++)
938  {
939  newkv[kind+i] = ub;
940  }
941  }
942  }
943  newkv.GetElements();
944 
945  swap(newpatch);
946 }
947 
948 void NURBSPatch::FlipDirection(int dir)
949 {
950  int size = SetLoopDirection(dir);
951 
952  for (int id = 0; id < nd/2; id++)
953  for (int i = 0; i < size; i++)
954  {
955  Swap<double>((*this)(id,i), (*this)(nd-1-id,i));
956  }
957  kv[dir]->Flip();
958 }
959 
960 void NURBSPatch::SwapDirections(int dir1, int dir2)
961 {
962  if (abs(dir1-dir2) == 2)
963  {
964  mfem_error("NURBSPatch::SwapDirections :"
965  " directions 0 and 2 are not supported!");
966  }
967 
968  Array<const KnotVector *> nkv(kv);
969 
970  Swap<const KnotVector *>(nkv[dir1], nkv[dir2]);
971  NURBSPatch *newpatch = new NURBSPatch(nkv, Dim);
972 
973  int size = SetLoopDirection(dir1);
974  newpatch->SetLoopDirection(dir2);
975 
976  for (int id = 0; id < nd; id++)
977  for (int i = 0; i < size; i++)
978  {
979  (*newpatch)(id,i) = (*this)(id,i);
980  }
981 
982  swap(newpatch);
983 }
984 
985 void NURBSPatch::Get3DRotationMatrix(double n[], double angle, double r,
986  DenseMatrix &T)
987 {
988  double c, s, c1;
989  double l2 = n[0]*n[0] + n[1]*n[1] + n[2]*n[2];
990  double l = sqrt(l2);
991 
992  if (fabs(angle) == M_PI_2)
993  {
994  s = r*copysign(1., angle);
995  c = 0.;
996  c1 = -1.;
997  }
998  else if (fabs(angle) == M_PI)
999  {
1000  s = 0.;
1001  c = -r;
1002  c1 = c - 1.;
1003  }
1004  else
1005  {
1006  s = r*sin(angle);
1007  c = r*cos(angle);
1008  c1 = c - 1.;
1009  }
1010 
1011  T.SetSize(3);
1012 
1013  T(0,0) = (n[0]*n[0] + (n[1]*n[1] + n[2]*n[2])*c)/l2;
1014  T(0,1) = -(n[0]*n[1]*c1)/l2 - (n[2]*s)/l;
1015  T(0,2) = -(n[0]*n[2]*c1)/l2 + (n[1]*s)/l;
1016  T(1,0) = -(n[0]*n[1]*c1)/l2 + (n[2]*s)/l;
1017  T(1,1) = (n[1]*n[1] + (n[0]*n[0] + n[2]*n[2])*c)/l2;
1018  T(1,2) = -(n[1]*n[2]*c1)/l2 - (n[0]*s)/l;
1019  T(2,0) = -(n[0]*n[2]*c1)/l2 - (n[1]*s)/l;
1020  T(2,1) = -(n[1]*n[2]*c1)/l2 + (n[0]*s)/l;
1021  T(2,2) = (n[2]*n[2] + (n[0]*n[0] + n[1]*n[1])*c)/l2;
1022 }
1023 
1024 void NURBSPatch::Rotate3D(double n[], double angle)
1025 {
1026  if (Dim != 4)
1027  {
1028  mfem_error("NURBSPatch::Rotate3D : not a NURBSPatch in 3D!");
1029  }
1030 
1031  DenseMatrix T(3);
1032  Vector x(3), y(NULL, 3);
1033 
1034  Get3DRotationMatrix(n, angle, 1., T);
1035 
1036  int size = 1;
1037  for (int i = 0; i < kv.Size(); i++)
1038  {
1039  size *= kv[i]->GetNCP();
1040  }
1041 
1042  for (int i = 0; i < size; i++)
1043  {
1044  y.SetData(data + i*Dim);
1045  x = y;
1046  T.Mult(x, y);
1047  }
1048 }
1049 
1050 int NURBSPatch::MakeUniformDegree(int degree)
1051 {
1052  int maxd = degree;
1053 
1054  if (maxd == -1)
1055  {
1056  for (int dir = 0; dir < kv.Size(); dir++)
1057  {
1058  maxd = std::max(maxd, kv[dir]->GetOrder());
1059  }
1060  }
1061 
1062  for (int dir = 0; dir < kv.Size(); dir++)
1063  {
1064  if (maxd > kv[dir]->GetOrder())
1065  {
1066  DegreeElevate(dir, maxd - kv[dir]->GetOrder());
1067  }
1068  }
1069 
1070  return maxd;
1071 }
1072 
1073 NURBSPatch *Interpolate(NURBSPatch &p1, NURBSPatch &p2)
1074 {
1075  if (p1.kv.Size() != p2.kv.Size() || p1.Dim != p2.Dim)
1076  {
1077  mfem_error("Interpolate(NURBSPatch &, NURBSPatch &)");
1078  }
1079 
1080  int size = 1, dim = p1.Dim;
1081  Array<const KnotVector *> kv(p1.kv.Size() + 1);
1082 
1083  for (int i = 0; i < p1.kv.Size(); i++)
1084  {
1085  if (p1.kv[i]->GetOrder() < p2.kv[i]->GetOrder())
1086  {
1087  p1.KnotInsert(i, *p2.kv[i]);
1088  p2.KnotInsert(i, *p1.kv[i]);
1089  }
1090  else
1091  {
1092  p2.KnotInsert(i, *p1.kv[i]);
1093  p1.KnotInsert(i, *p2.kv[i]);
1094  }
1095  kv[i] = p1.kv[i];
1096  size *= kv[i]->GetNCP();
1097  }
1098 
1099  KnotVector &nkv = *(new KnotVector(1, 2));
1100  nkv[0] = nkv[1] = 0.0;
1101  nkv[2] = nkv[3] = 1.0;
1102  nkv.GetElements();
1103  kv.Last() = &nkv;
1104 
1105  NURBSPatch *patch = new NURBSPatch(kv, dim);
1106  delete kv.Last();
1107 
1108  for (int i = 0; i < size; i++)
1109  {
1110  for (int d = 0; d < dim; d++)
1111  {
1112  patch->data[i*dim+d] = p1.data[i*dim+d];
1113  patch->data[(i+size)*dim+d] = p2.data[i*dim+d];
1114  }
1115  }
1116 
1117  return patch;
1118 }
1119 
1120 NURBSPatch *Revolve3D(NURBSPatch &patch, double n[], double ang, int times)
1121 {
1122  if (patch.Dim != 4)
1123  {
1124  mfem_error("Revolve3D(NURBSPatch &, double [], double)");
1125  }
1126 
1127  int size = 1, ns;
1128  Array<const KnotVector *> nkv(patch.kv.Size() + 1);
1129 
1130  for (int i = 0; i < patch.kv.Size(); i++)
1131  {
1132  nkv[i] = patch.kv[i];
1133  size *= nkv[i]->GetNCP();
1134  }
1135  ns = 2*times + 1;
1136  KnotVector &lkv = *(new KnotVector(2, ns));
1137  nkv.Last() = &lkv;
1138  lkv[0] = lkv[1] = lkv[2] = 0.0;
1139  for (int i = 1; i < times; i++)
1140  {
1141  lkv[2*i+1] = lkv[2*i+2] = i;
1142  }
1143  lkv[ns] = lkv[ns+1] = lkv[ns+2] = times;
1144  lkv.GetElements();
1145  NURBSPatch *newpatch = new NURBSPatch(nkv, 4);
1146  delete nkv.Last();
1147 
1148  DenseMatrix T(3), T2(3);
1149  Vector u(NULL, 3), v(NULL, 3);
1150 
1151  NURBSPatch::Get3DRotationMatrix(n, ang, 1., T);
1152  double c = cos(ang/2);
1153  NURBSPatch::Get3DRotationMatrix(n, ang/2, 1./c, T2);
1154  T2 *= c;
1155 
1156  double *op = patch.data, *np;
1157  for (int i = 0; i < size; i++)
1158  {
1159  np = newpatch->data + 4*i;
1160  for (int j = 0; j < 4; j++)
1161  {
1162  np[j] = op[j];
1163  }
1164  for (int j = 0; j < times; j++)
1165  {
1166  u.SetData(np);
1167  v.SetData(np += 4*size);
1168  T2.Mult(u, v);
1169  v[3] = c*u[3];
1170  v.SetData(np += 4*size);
1171  T.Mult(u, v);
1172  v[3] = u[3];
1173  }
1174  op += 4;
1175  }
1176 
1177  return newpatch;
1178 }
1179 
1180 
1181 NURBSExtension::NURBSExtension(const NURBSExtension &orig)
1182  : mOrder(orig.mOrder), mOrders(orig.mOrders),
1183  NumOfKnotVectors(orig.NumOfKnotVectors),
1184  NumOfVertices(orig.NumOfVertices),
1185  NumOfElements(orig.NumOfElements),
1186  NumOfBdrElements(orig.NumOfBdrElements),
1187  NumOfDofs(orig.NumOfDofs),
1188  NumOfActiveVertices(orig.NumOfActiveVertices),
1189  NumOfActiveElems(orig.NumOfActiveElems),
1190  NumOfActiveBdrElems(orig.NumOfActiveBdrElems),
1191  NumOfActiveDofs(orig.NumOfActiveDofs),
1192  activeVert(orig.activeVert),
1193  activeElem(orig.activeElem),
1194  activeBdrElem(orig.activeBdrElem),
1195  activeDof(orig.activeDof),
1196  patchTopo(new Mesh(*orig.patchTopo)),
1197  own_topo(true),
1198  edge_to_knot(orig.edge_to_knot),
1199  knotVectors(orig.knotVectors.Size()), // knotVectors are copied in the body
1200  weights(orig.weights),
1201  d_to_d(orig.d_to_d),
1202  master(orig.master),
1203  slave(orig.slave),
1204  v_meshOffsets(orig.v_meshOffsets),
1205  e_meshOffsets(orig.e_meshOffsets),
1206  f_meshOffsets(orig.f_meshOffsets),
1207  p_meshOffsets(orig.p_meshOffsets),
1208  v_spaceOffsets(orig.v_spaceOffsets),
1209  e_spaceOffsets(orig.e_spaceOffsets),
1210  f_spaceOffsets(orig.f_spaceOffsets),
1211  p_spaceOffsets(orig.p_spaceOffsets),
1212  el_dof(orig.el_dof ? new Table(*orig.el_dof) : NULL),
1213  bel_dof(orig.bel_dof ? new Table(*orig.bel_dof) : NULL),
1214  el_to_patch(orig.el_to_patch),
1215  bel_to_patch(orig.bel_to_patch),
1216  el_to_IJK(orig.el_to_IJK),
1217  bel_to_IJK(orig.bel_to_IJK),
1218  patches(orig.patches.Size()) // patches are copied in the body
1219 {
1220  // Copy the knot vectors:
1221  for (int i = 0; i < knotVectors.Size(); i++)
1222  {
1223  knotVectors[i] = new KnotVector(*orig.knotVectors[i]);
1224  }
1225 
1226  // Copy the patches:
1227  for (int p = 0; p < patches.Size(); p++)
1228  {
1229  patches[p] = new NURBSPatch(*orig.patches[p]);
1230  }
1231 }
1232 
1233 NURBSExtension::NURBSExtension(std::istream &input)
1234 {
1235  // Read topology
1236  patchTopo = new Mesh;
1237  patchTopo->LoadPatchTopo(input, edge_to_knot);
1238  own_topo = 1;
1239 
1240  CheckPatches();
1241  // CheckBdrPatches();
1242 
1243  skip_comment_lines(input, '#');
1244 
1245  // Read knotvectors or patches
1246  string ident;
1247  input >> ws >> ident; // 'knotvectors' or 'patches'
1248  if (ident == "knotvectors")
1249  {
1250  input >> NumOfKnotVectors;
1251  knotVectors.SetSize(NumOfKnotVectors);
1252  for (int i = 0; i < NumOfKnotVectors; i++)
1253  {
1254  knotVectors[i] = new KnotVector(input);
1255  }
1256  }
1257  else if (ident == "patches")
1258  {
1259  patches.SetSize(GetNP());
1260  for (int p = 0; p < patches.Size(); p++)
1261  {
1262  skip_comment_lines(input, '#');
1263  patches[p] = new NURBSPatch(input);
1264  }
1265 
1266  NumOfKnotVectors = 0;
1267  for (int i = 0; i < patchTopo->GetNEdges(); i++)
1268  if (NumOfKnotVectors < KnotInd(i))
1269  {
1270  NumOfKnotVectors = KnotInd(i);
1271  }
1272  NumOfKnotVectors++;
1273  knotVectors.SetSize(NumOfKnotVectors);
1274  knotVectors = NULL;
1275 
1276  Array<int> edges, oedge;
1277  for (int p = 0; p < patches.Size(); p++)
1278  {
1279  patchTopo->GetElementEdges(p, edges, oedge);
1280  if (Dimension() == 2)
1281  {
1282  if (knotVectors[KnotInd(edges[0])] == NULL)
1283  {
1284  knotVectors[KnotInd(edges[0])] =
1285  new KnotVector(*patches[p]->GetKV(0));
1286  }
1287  if (knotVectors[KnotInd(edges[1])] == NULL)
1288  {
1289  knotVectors[KnotInd(edges[1])] =
1290  new KnotVector(*patches[p]->GetKV(1));
1291  }
1292  }
1293  else
1294  {
1295  if (knotVectors[KnotInd(edges[0])] == NULL)
1296  {
1297  knotVectors[KnotInd(edges[0])] =
1298  new KnotVector(*patches[p]->GetKV(0));
1299  }
1300  if (knotVectors[KnotInd(edges[3])] == NULL)
1301  {
1302  knotVectors[KnotInd(edges[3])] =
1303  new KnotVector(*patches[p]->GetKV(1));
1304  }
1305  if (knotVectors[KnotInd(edges[8])] == NULL)
1306  {
1307  knotVectors[KnotInd(edges[8])] =
1308  new KnotVector(*patches[p]->GetKV(2));
1309  }
1310  }
1311  }
1312  }
1313  else
1314  {
1315  MFEM_ABORT("invalid section: " << ident);
1316  }
1317 
1318  SetOrdersFromKnotVectors();
1319 
1320  GenerateOffsets();
1321  CountElements();
1322  CountBdrElements();
1323  // NumOfVertices, NumOfElements, NumOfBdrElements, NumOfDofs
1324 
1325  skip_comment_lines(input, '#');
1326 
1327  // Check for a list of mesh elements
1328  if (patches.Size() == 0)
1329  {
1330  input >> ws >> ident;
1331  }
1332  if (patches.Size() == 0 && ident == "mesh_elements")
1333  {
1334  input >> NumOfActiveElems;
1335  activeElem.SetSize(GetGNE());
1336  activeElem = false;
1337  int glob_elem;
1338  for (int i = 0; i < NumOfActiveElems; i++)
1339  {
1340  input >> glob_elem;
1341  activeElem[glob_elem] = true;
1342  }
1343 
1344  skip_comment_lines(input, '#');
1345  input >> ws >> ident;
1346  }
1347  else
1348  {
1349  NumOfActiveElems = NumOfElements;
1350  activeElem.SetSize(NumOfElements);
1351  activeElem = true;
1352  }
1353 
1354  GenerateActiveVertices();
1355  InitDofMap();
1356  GenerateElementDofTable();
1357  GenerateActiveBdrElems();
1358  GenerateBdrElementDofTable();
1359 
1360  // periodic
1361  if (ident == "periodic")
1362  {
1363  master.Load(input);
1364  slave.Load(input);
1365 
1366  skip_comment_lines(input, '#');
1367  input >> ws >> ident;
1368  }
1369 
1370  if (patches.Size() == 0)
1371  {
1372  // weights
1373  if (ident == "weights")
1374  {
1375  weights.Load(input, GetNDof());
1376  }
1377  else // e.g. ident = "unitweights" or "autoweights"
1378  {
1379  weights.SetSize(GetNDof());
1380  weights = 1.0;
1381  }
1382  }
1383 
1384  // periodic
1385  ConnectBoundaries();
1386 }
1387 
1388 NURBSExtension::NURBSExtension(NURBSExtension *parent, int newOrder)
1389 {
1390  patchTopo = parent->patchTopo;
1391  own_topo = 0;
1392 
1393  parent->edge_to_knot.Copy(edge_to_knot);
1394 
1395  NumOfKnotVectors = parent->GetNKV();
1396  knotVectors.SetSize(NumOfKnotVectors);
1397  const Array<int> &pOrders = parent->GetOrders();
1398  for (int i = 0; i < NumOfKnotVectors; i++)
1399  {
1400  if (newOrder > pOrders[i])
1401  {
1402  knotVectors[i] =
1403  parent->GetKnotVector(i)->DegreeElevate(newOrder - pOrders[i]);
1404  }
1405  else
1406  {
1407  knotVectors[i] = new KnotVector(*parent->GetKnotVector(i));
1408  }
1409  }
1410 
1411  // copy some data from parent
1412  NumOfElements = parent->NumOfElements;
1413  NumOfBdrElements = parent->NumOfBdrElements;
1414 
1415  SetOrdersFromKnotVectors();
1416 
1417  GenerateOffsets(); // dof offsets will be different from parent
1418 
1419  NumOfActiveVertices = parent->NumOfActiveVertices;
1420  NumOfActiveElems = parent->NumOfActiveElems;
1421  NumOfActiveBdrElems = parent->NumOfActiveBdrElems;
1422  parent->activeVert.Copy(activeVert);
1423  InitDofMap();
1424  parent->activeElem.Copy(activeElem);
1425  parent->activeBdrElem.Copy(activeBdrElem);
1426 
1427  GenerateElementDofTable();
1428  GenerateBdrElementDofTable();
1429 
1430  weights.SetSize(GetNDof());
1431  weights = 1.0;
1432 
1433  // periodic
1434  parent->master.Copy(master);
1435  parent->slave.Copy(slave);
1436  ConnectBoundaries();
1437 }
1438 
1439 NURBSExtension::NURBSExtension(NURBSExtension *parent,
1440  const Array<int> &newOrders)
1441 {
1442  newOrders.Copy(mOrders);
1443  SetOrderFromOrders();
1444 
1445  patchTopo = parent->patchTopo;
1446  own_topo = 0;
1447 
1448  parent->edge_to_knot.Copy(edge_to_knot);
1449 
1450  NumOfKnotVectors = parent->GetNKV();
1451  MFEM_VERIFY(mOrders.Size() == NumOfKnotVectors, "invalid newOrders array");
1452  knotVectors.SetSize(NumOfKnotVectors);
1453  const Array<int> &pOrders = parent->GetOrders();
1454 
1455  for (int i = 0; i < NumOfKnotVectors; i++)
1456  {
1457  if (mOrders[i] > pOrders[i])
1458  {
1459  knotVectors[i] =
1460  parent->GetKnotVector(i)->DegreeElevate(mOrders[i] - pOrders[i]);
1461  }
1462  else
1463  {
1464  knotVectors[i] = new KnotVector(*parent->GetKnotVector(i));
1465  }
1466  }
1467 
1468  // copy some data from parent
1469  NumOfElements = parent->NumOfElements;
1470  NumOfBdrElements = parent->NumOfBdrElements;
1471 
1472  GenerateOffsets(); // dof offsets will be different from parent
1473 
1474  NumOfActiveVertices = parent->NumOfActiveVertices;
1475  NumOfActiveElems = parent->NumOfActiveElems;
1476  NumOfActiveBdrElems = parent->NumOfActiveBdrElems;
1477  parent->activeVert.Copy(activeVert);
1478  InitDofMap();
1479  parent->activeElem.Copy(activeElem);
1480  parent->activeBdrElem.Copy(activeBdrElem);
1481 
1482  GenerateElementDofTable();
1483  GenerateBdrElementDofTable();
1484 
1485  weights.SetSize(GetNDof());
1486  weights = 1.0;
1487 
1488  parent->master.Copy(master);
1489  parent->slave.Copy(slave);
1490  ConnectBoundaries();
1491 }
1492 
1493 NURBSExtension::NURBSExtension(Mesh *mesh_array[], int num_pieces)
1494 {
1495  NURBSExtension *parent = mesh_array[0]->NURBSext;
1496 
1497  if (!parent->own_topo)
1498  {
1499  mfem_error("NURBSExtension::NURBSExtension :\n"
1500  " parent does not own the patch topology!");
1501  }
1502  patchTopo = parent->patchTopo;
1503  own_topo = 1;
1504  parent->own_topo = 0;
1505 
1506  parent->edge_to_knot.Copy(edge_to_knot);
1507 
1508  parent->GetOrders().Copy(mOrders);
1509  mOrder = parent->GetOrder();
1510 
1511  NumOfKnotVectors = parent->GetNKV();
1512  knotVectors.SetSize(NumOfKnotVectors);
1513  for (int i = 0; i < NumOfKnotVectors; i++)
1514  {
1515  knotVectors[i] = new KnotVector(*parent->GetKnotVector(i));
1516  }
1517 
1518  GenerateOffsets();
1519  CountElements();
1520  CountBdrElements();
1521 
1522  // assuming the meshes define a partitioning of all the elements
1523  NumOfActiveElems = NumOfElements;
1524  activeElem.SetSize(NumOfElements);
1525  activeElem = true;
1526 
1527  GenerateActiveVertices();
1528  InitDofMap();
1529  GenerateElementDofTable();
1530  GenerateActiveBdrElems();
1531  GenerateBdrElementDofTable();
1532 
1533  weights.SetSize(GetNDof());
1534  MergeWeights(mesh_array, num_pieces);
1535 }
1536 
1537 NURBSExtension::~NURBSExtension()
1538 {
1539  if (patches.Size() == 0)
1540  {
1541  delete bel_dof;
1542  delete el_dof;
1543  }
1544 
1545  for (int i = 0; i < knotVectors.Size(); i++)
1546  {
1547  delete knotVectors[i];
1548  }
1549 
1550  for (int i = 0; i < patches.Size(); i++)
1551  {
1552  delete patches[i];
1553  }
1554 
1555  if (own_topo)
1556  {
1557  delete patchTopo;
1558  }
1559 }
1560 
1561 void NURBSExtension::Print(std::ostream &out) const
1562 {
1563  patchTopo->PrintTopo(out, edge_to_knot);
1564  if (patches.Size() == 0)
1565  {
1566  out << "\nknotvectors\n" << NumOfKnotVectors << '\n';
1567  for (int i = 0; i < NumOfKnotVectors; i++)
1568  {
1569  knotVectors[i]->Print(out);
1570  }
1571 
1572  if (NumOfActiveElems < NumOfElements)
1573  {
1574  out << "\nmesh_elements\n" << NumOfActiveElems << '\n';
1575  for (int i = 0; i < NumOfElements; i++)
1576  if (activeElem[i])
1577  {
1578  out << i << '\n';
1579  }
1580  }
1581 
1582  out << "\nweights\n";
1583  weights.Print(out, 1);
1584  }
1585  else
1586  {
1587  out << "\npatches\n";
1588  for (int p = 0; p < patches.Size(); p++)
1589  {
1590  out << "\n# patch " << p << "\n\n";
1591  patches[p]->Print(out);
1592  }
1593  }
1594 }
1595 
1596 void NURBSExtension::PrintCharacteristics(std::ostream &out) const
1597 {
1598  out <<
1599  "NURBS Mesh entity sizes:\n"
1600  "Dimension = " << Dimension() << "\n"
1601  "Unique Orders = ";
1602  Array<int> unique_orders(mOrders);
1603  unique_orders.Sort();
1604  unique_orders.Unique();
1605  unique_orders.Print(out, unique_orders.Size());
1606  out <<
1607  "NumOfKnotVectors = " << GetNKV() << "\n"
1608  "NumOfPatches = " << GetNP() << "\n"
1609  "NumOfBdrPatches = " << GetNBP() << "\n"
1610  "NumOfVertices = " << GetGNV() << "\n"
1611  "NumOfElements = " << GetGNE() << "\n"
1612  "NumOfBdrElements = " << GetGNBE() << "\n"
1613  "NumOfDofs = " << GetNTotalDof() << "\n"
1614  "NumOfActiveVertices = " << GetNV() << "\n"
1615  "NumOfActiveElems = " << GetNE() << "\n"
1616  "NumOfActiveBdrElems = " << GetNBE() << "\n"
1617  "NumOfActiveDofs = " << GetNDof() << '\n';
1618  for (int i = 0; i < NumOfKnotVectors; i++)
1619  {
1620  out << ' ' << i + 1 << ") ";
1621  knotVectors[i]->Print(out);
1622  }
1623  out << endl;
1624 }
1625 
1626 void NURBSExtension::InitDofMap()
1627 {
1628  master.SetSize(0);
1629  slave.SetSize(0);
1630  d_to_d.SetSize(0);
1631 }
1632 
1633 void NURBSExtension::ConnectBoundaries(Array<int> &bnds0, Array<int> &bnds1)
1634 {
1635  bnds0.Copy(master);
1636  bnds1.Copy(slave);
1637  ConnectBoundaries();
1638 }
1639 
1640 void NURBSExtension::ConnectBoundaries()
1641 {
1642  if (master.Size() != slave.Size())
1643  {
1644  mfem_error("NURBSExtension::ConnectBoundaries() boundary lists not of equal size");
1645  }
1646  if (master.Size() == 0 ) { return; }
1647 
1648  // Connect
1649  d_to_d.SetSize(NumOfDofs);
1650  for (int i = 0; i < NumOfDofs; i++) { d_to_d[i] = i; }
1651 
1652  for (int i = 0; i < master.Size(); i++)
1653  {
1654  if (Dimension() == 2)
1655  {
1656  ConnectBoundaries2D(master[i], slave[i]);
1657  }
1658  else
1659  {
1660  ConnectBoundaries3D(master[i], slave[i]);
1661  }
1662  }
1663 
1664  // Finalize
1665  GenerateElementDofTable();
1666  GenerateBdrElementDofTable();
1667 }
1668 
1669 void NURBSExtension::ConnectBoundaries2D(int bnd0, int bnd1)
1670 {
1671  int idx0 = -1, idx1 = -1;
1672  for (int b = 0; b < GetNBP(); b++)
1673  {
1674  if (bnd0 == patchTopo->GetBdrAttribute(b)) { idx0 = b; }
1675  if (bnd1 == patchTopo->GetBdrAttribute(b)) { idx1 = b; }
1676  }
1677  MFEM_VERIFY(idx0 != -1,"Bdr 0 not found");
1678  MFEM_VERIFY(idx1 != -1,"Bdr 1 not found");
1679 
1680  NURBSPatchMap p2g0(this);
1681  NURBSPatchMap p2g1(this);
1682 
1683  int okv0[1],okv1[1];
1684  const KnotVector *kv0[1],*kv1[1];
1685 
1686  p2g0.SetBdrPatchDofMap(idx0, kv0, okv0);
1687  p2g1.SetBdrPatchDofMap(idx1, kv1, okv1);
1688 
1689  int nx = p2g0.nx();
1690  int nks0 = kv0[0]->GetNKS();
1691 
1692 #ifdef MFEM_DEBUG
1693  bool compatible = true;
1694  if (p2g0.nx() != p2g1.nx()) { compatible = false; }
1695  if (kv0[0]->GetNKS() != kv1[0]->GetNKS()) { compatible = false; }
1696  if (kv0[0]->GetOrder() != kv1[0]->GetOrder()) { compatible = false; }
1697 
1698  if (!compatible)
1699  {
1700  mfem::out<<p2g0.nx()<<" "<<p2g1.nx()<<endl;
1701  mfem::out<<kv0[0]->GetNKS()<<" "<<kv1[0]->GetNKS()<<endl;
1702  mfem::out<<kv0[0]->GetOrder()<<" "<<kv1[0]->GetOrder()<<endl;
1703  mfem_error("NURBS boundaries not compatible");
1704  }
1705 #endif
1706 
1707  for (int i = 0; i < nks0; i++)
1708  {
1709  if (kv0[0]->isElement(i))
1710  {
1711  if (!kv1[0]->isElement(i)) { mfem_error("isElement does not match"); }
1712  for (int ii = 0; ii <= kv0[0]->GetOrder(); ii++)
1713  {
1714  int ii0 = (okv0[0] >= 0) ? (i+ii) : (nx-i-ii);
1715  int ii1 = (okv1[0] >= 0) ? (i+ii) : (nx-i-ii);
1716 
1717  d_to_d[p2g0(ii0)] = d_to_d[p2g1(ii1)];
1718  }
1719 
1720  }
1721  }
1722 
1723  // Clean d_to_d
1724  Array<int> tmp(d_to_d.Size()+1);
1725  tmp = 0;
1726 
1727  for (int i = 0; i < d_to_d.Size(); i++)
1728  {
1729  tmp[d_to_d[i]] = 1;
1730  }
1731 
1732  int cnt = 0;
1733  for (int i = 0; i < tmp.Size(); i++)
1734  {
1735  if (tmp[i] == 1) { tmp[i] = cnt++; }
1736  }
1737  NumOfDofs = cnt;
1738 
1739  for (int i = 0; i < d_to_d.Size(); i++)
1740  {
1741  d_to_d[i] = tmp[d_to_d[i]];
1742  }
1743 }
1744 
1745 void NURBSExtension::ConnectBoundaries3D(int bnd0, int bnd1)
1746 {
1747  int idx0 = -1, idx1 = -1;
1748  for (int b = 0; b < GetNBP(); b++)
1749  {
1750  if (bnd0 == patchTopo->GetBdrAttribute(b)) { idx0 = b; }
1751  if (bnd1 == patchTopo->GetBdrAttribute(b)) { idx1 = b; }
1752  }
1753  MFEM_VERIFY(idx0 != -1,"Bdr 0 not found");
1754  MFEM_VERIFY(idx1 != -1,"Bdr 1 not found");
1755 
1756  NURBSPatchMap p2g0(this);
1757  NURBSPatchMap p2g1(this);
1758 
1759  int okv0[2],okv1[2];
1760  const KnotVector *kv0[2],*kv1[2];
1761 
1762  p2g0.SetBdrPatchDofMap(idx0, kv0, okv0);
1763  p2g1.SetBdrPatchDofMap(idx1, kv1, okv1);
1764 
1765  int nx = p2g0.nx();
1766  int ny = p2g0.ny();
1767 
1768  int nks0 = kv0[0]->GetNKS();
1769  int nks1 = kv0[1]->GetNKS();
1770 
1771 #ifdef MFEM_DEBUG
1772  bool compatible = true;
1773  if (p2g0.nx() != p2g1.nx()) { compatible = false; }
1774  if (p2g0.ny() != p2g1.ny()) { compatible = false; }
1775 
1776  if (kv0[0]->GetNKS() != kv1[0]->GetNKS()) { compatible = false; }
1777  if (kv0[1]->GetNKS() != kv1[0]->GetNKS()) { compatible = false; }
1778 
1779  if (kv0[0]->GetOrder() != kv1[0]->GetOrder()) { compatible = false; }
1780  if (kv0[1]->GetOrder() != kv1[1]->GetOrder()) { compatible = false; }
1781 
1782  if (!compatible)
1783  {
1784  mfem::out<<p2g0.nx()<<" "<<p2g1.nx()<<endl;
1785  mfem::out<<p2g0.ny()<<" "<<p2g1.ny()<<endl;
1786 
1787  mfem::out<<kv0[0]->GetNKS()<<" "<<kv1[0]->GetNKS()<<endl;
1788  mfem::out<<kv0[1]->GetNKS()<<" "<<kv1[0]->GetNKS()<<endl;
1789 
1790  mfem::out<<kv0[0]->GetOrder()<<" "<<kv1[0]->GetOrder()<<endl;
1791  mfem::out<<kv0[1]->GetOrder()<<" "<<kv1[1]->GetOrder()<<endl;
1792  mfem_error("NURBS boundaries not compatible");
1793  }
1794 #endif
1795 
1796  for (int j = 0; j < nks1; j++)
1797  {
1798  if (kv0[1]->isElement(j))
1799  {
1800  if (!kv1[1]->isElement(j)) { mfem_error("isElement does not match #1"); }
1801  for (int i = 0; i < nks0; i++)
1802  {
1803  if (kv0[0]->isElement(i))
1804  {
1805  if (!kv1[0]->isElement(i)) { mfem_error("isElement does not match #0"); }
1806  for (int jj = 0; jj <= kv0[1]->GetOrder(); jj++)
1807  {
1808  int jj0 = (okv0[1] >= 0) ? (j+jj) : (ny-j-jj);
1809  int jj1 = (okv1[1] >= 0) ? (j+jj) : (ny-j-jj);
1810 
1811  for (int ii = 0; ii <= kv0[0]->GetOrder(); ii++)
1812  {
1813  int ii0 = (okv0[0] >= 0) ? (i+ii) : (nx-i-ii);
1814  int ii1 = (okv1[0] >= 0) ? (i+ii) : (nx-i-ii);
1815 
1816  d_to_d[p2g0(ii0,jj0)] = d_to_d[p2g1(ii1,jj1)];
1817  }
1818  }
1819  }
1820  }
1821  }
1822  }
1823 
1824  // Clean d_to_d
1825  Array<int> tmp(d_to_d.Size()+1);
1826  tmp = 0;
1827 
1828  for (int i = 0; i < d_to_d.Size(); i++)
1829  {
1830  tmp[d_to_d[i]] = 1;
1831  }
1832 
1833  int cnt = 0;
1834  for (int i = 0; i < tmp.Size(); i++)
1835  {
1836  if (tmp[i] == 1) { tmp[i] = cnt++; }
1837  }
1838  NumOfDofs = cnt;
1839 
1840  for (int i = 0; i < d_to_d.Size(); i++)
1841  {
1842  d_to_d[i] = tmp[d_to_d[i]];
1843  }
1844 }
1845 
1847 {
1848  int vert[8], nv, g_el, nx, ny, nz, dim = Dimension();
1849 
1850  NURBSPatchMap p2g(this);
1851  const KnotVector *kv[3];
1852 
1853  g_el = 0;
1854  activeVert.SetSize(GetGNV());
1855  activeVert = -1;
1856  for (int p = 0; p < GetNP(); p++)
1857  {
1858  p2g.SetPatchVertexMap(p, kv);
1859 
1860  nx = p2g.nx();
1861  ny = p2g.ny();
1862  nz = (dim == 3) ? p2g.nz() : 1;
1863 
1864  for (int k = 0; k < nz; k++)
1865  {
1866  for (int j = 0; j < ny; j++)
1867  {
1868  for (int i = 0; i < nx; i++)
1869  {
1870  if (activeElem[g_el])
1871  {
1872  if (dim == 2)
1873  {
1874  vert[0] = p2g(i, j );
1875  vert[1] = p2g(i+1,j );
1876  vert[2] = p2g(i+1,j+1);
1877  vert[3] = p2g(i, j+1);
1878  nv = 4;
1879  }
1880  else
1881  {
1882  vert[0] = p2g(i, j, k);
1883  vert[1] = p2g(i+1,j, k);
1884  vert[2] = p2g(i+1,j+1,k);
1885  vert[3] = p2g(i, j+1,k);
1886 
1887  vert[4] = p2g(i, j, k+1);
1888  vert[5] = p2g(i+1,j, k+1);
1889  vert[6] = p2g(i+1,j+1,k+1);
1890  vert[7] = p2g(i, j+1,k+1);
1891  nv = 8;
1892  }
1893 
1894  for (int v = 0; v < nv; v++)
1895  {
1896  activeVert[vert[v]] = 1;
1897  }
1898  }
1899  g_el++;
1900  }
1901  }
1902  }
1903  }
1904 
1905  NumOfActiveVertices = 0;
1906  for (int i = 0; i < GetGNV(); i++)
1907  if (activeVert[i] == 1)
1908  {
1909  activeVert[i] = NumOfActiveVertices++;
1910  }
1911 }
1912 
1914 {
1915  int dim = Dimension();
1916  Array<KnotVector *> kv(dim);
1917 
1918  activeBdrElem.SetSize(GetGNBE());
1919  if (GetGNE() == GetNE())
1920  {
1921  activeBdrElem = true;
1922  NumOfActiveBdrElems = GetGNBE();
1923  return;
1924  }
1925  activeBdrElem = false;
1926  NumOfActiveBdrElems = 0;
1927  // the mesh will generate the actual boundary including boundary
1928  // elements that are not on boundary patches. we use this for
1929  // visualization of processor boundaries
1930 
1931  // TODO: generate actual boundary?
1932 }
1933 
1934 
1935 void NURBSExtension::MergeWeights(Mesh *mesh_array[], int num_pieces)
1936 {
1937  Array<int> lelem_elem;
1938 
1939  for (int i = 0; i < num_pieces; i++)
1940  {
1941  NURBSExtension *lext = mesh_array[i]->NURBSext;
1942 
1943  lext->GetElementLocalToGlobal(lelem_elem);
1944 
1945  for (int lel = 0; lel < lext->GetNE(); lel++)
1946  {
1947  int gel = lelem_elem[lel];
1948 
1949  int nd = el_dof->RowSize(gel);
1950  int *gdofs = el_dof->GetRow(gel);
1951  int *ldofs = lext->el_dof->GetRow(lel);
1952  for (int j = 0; j < nd; j++)
1953  {
1954  weights(gdofs[j]) = lext->weights(ldofs[j]);
1955  }
1956  }
1957  }
1958 }
1959 
1961  GridFunction *gf_array[], int num_pieces, GridFunction &merged)
1962 {
1963  FiniteElementSpace *gfes = merged.FESpace();
1964  Array<int> lelem_elem, dofs;
1965  Vector lvec;
1966 
1967  for (int i = 0; i < num_pieces; i++)
1968  {
1969  FiniteElementSpace *lfes = gf_array[i]->FESpace();
1970  NURBSExtension *lext = lfes->GetMesh()->NURBSext;
1971 
1972  lext->GetElementLocalToGlobal(lelem_elem);
1973 
1974  for (int lel = 0; lel < lext->GetNE(); lel++)
1975  {
1976  lfes->GetElementVDofs(lel, dofs);
1977  gf_array[i]->GetSubVector(dofs, lvec);
1978 
1979  gfes->GetElementVDofs(lelem_elem[lel], dofs);
1980  merged.SetSubVector(dofs, lvec);
1981  }
1982  }
1983 }
1984 
1986 {
1987  Array<int> edges;
1988  Array<int> oedge;
1989 
1990  for (int p = 0; p < GetNP(); p++)
1991  {
1992  patchTopo->GetElementEdges(p, edges, oedge);
1993 
1994  for (int i = 0; i < edges.Size(); i++)
1995  {
1996  edges[i] = edge_to_knot[edges[i]];
1997  if (oedge[i] < 0)
1998  {
1999  edges[i] = -1 - edges[i];
2000  }
2001  }
2002 
2003  if ((Dimension() == 2 &&
2004  (edges[0] != -1 - edges[2] || edges[1] != -1 - edges[3])) ||
2005 
2006  (Dimension() == 3 &&
2007  (edges[0] != edges[2] || edges[0] != edges[4] ||
2008  edges[0] != edges[6] || edges[1] != edges[3] ||
2009  edges[1] != edges[5] || edges[1] != edges[7] ||
2010  edges[8] != edges[9] || edges[8] != edges[10] ||
2011  edges[8] != edges[11])))
2012  {
2013  mfem::err << "NURBSExtension::CheckPatch (patch = " << p
2014  << ")\n Inconsistent edge-to-knot mapping!\n";
2015  mfem_error();
2016  }
2017 
2018  if ((Dimension() == 2 &&
2019  (edges[0] < 0 || edges[1] < 0)) ||
2020 
2021  (Dimension() == 3 &&
2022  (edges[0] < 0 || edges[3] < 0 || edges[8] < 0)))
2023  {
2024  mfem::err << "NURBSExtension::CheckPatch (patch = " << p
2025  << ") : Bad orientation!\n";
2026  mfem_error();
2027  }
2028  }
2029 }
2030 
2032 {
2033  Array<int> edges;
2034  Array<int> oedge;
2035 
2036  for (int p = 0; p < GetNBP(); p++)
2037  {
2038  patchTopo->GetBdrElementEdges(p, edges, oedge);
2039 
2040  for (int i = 0; i < edges.Size(); i++)
2041  {
2042  edges[i] = edge_to_knot[edges[i]];
2043  if (oedge[i] < 0)
2044  {
2045  edges[i] = -1 - edges[i];
2046  }
2047  }
2048 
2049  if ((Dimension() == 2 && (edges[0] < 0)) ||
2050  (Dimension() == 3 && (edges[0] < 0 || edges[1] < 0)))
2051  {
2052  mfem::err << "NURBSExtension::CheckBdrPatch (boundary patch = "
2053  << p << ") : Bad orientation!\n";
2054  mfem_error();
2055  }
2056  }
2057 }
2058 
2060 {
2061  Array<int> edges, orient;
2062 
2063  kv.SetSize(Dimension());
2064  patchTopo->GetElementEdges(p, edges, orient);
2065 
2066  if (Dimension() == 2)
2067  {
2068  kv[0] = KnotVec(edges[0]);
2069  kv[1] = KnotVec(edges[1]);
2070  }
2071  else
2072  {
2073  kv[0] = KnotVec(edges[0]);
2074  kv[1] = KnotVec(edges[3]);
2075  kv[2] = KnotVec(edges[8]);
2076  }
2077 }
2078 
2080 const
2081 {
2082  Array<int> edges, orient;
2083 
2084  kv.SetSize(Dimension());
2085  patchTopo->GetElementEdges(p, edges, orient);
2086 
2087  if (Dimension() == 2)
2088  {
2089  kv[0] = KnotVec(edges[0]);
2090  kv[1] = KnotVec(edges[1]);
2091  }
2092  else
2093  {
2094  kv[0] = KnotVec(edges[0]);
2095  kv[1] = KnotVec(edges[3]);
2096  kv[2] = KnotVec(edges[8]);
2097  }
2098 }
2099 
2101 {
2102  Array<int> edges;
2103  Array<int> orient;
2104 
2105  kv.SetSize(Dimension() - 1);
2106  patchTopo->GetBdrElementEdges(p, edges, orient);
2107 
2108  if (Dimension() == 2)
2109  {
2110  kv[0] = KnotVec(edges[0]);
2111  }
2112  else
2113  {
2114  kv[0] = KnotVec(edges[0]);
2115  kv[1] = KnotVec(edges[1]);
2116  }
2117 }
2118 
2120  int p, Array<const KnotVector *> &kv) const
2121 {
2122  Array<int> edges;
2123  Array<int> orient;
2124 
2125  kv.SetSize(Dimension() - 1);
2126  patchTopo->GetBdrElementEdges(p, edges, orient);
2127 
2128  if (Dimension() == 2)
2129  {
2130  kv[0] = KnotVec(edges[0]);
2131  }
2132  else
2133  {
2134  kv[0] = KnotVec(edges[0]);
2135  kv[1] = KnotVec(edges[1]);
2136  }
2137 }
2138 
2140 {
2141  MFEM_VERIFY(mOrders.Size() > 0, "");
2142  mOrder = mOrders[0];
2143  for (int i = 1; i < mOrders.Size(); i++)
2144  {
2145  if (mOrders[i] != mOrder)
2146  {
2148  return;
2149  }
2150  }
2151 }
2152 
2154 {
2155  mOrders.SetSize(NumOfKnotVectors);
2156  for (int i = 0; i < NumOfKnotVectors; i++)
2157  {
2158  mOrders[i] = knotVectors[i]->GetOrder();
2159  }
2160  SetOrderFromOrders();
2161 }
2162 
2164 {
2165  int nv = patchTopo->GetNV();
2166  int ne = patchTopo->GetNEdges();
2167  int nf = patchTopo->GetNFaces();
2168  int np = patchTopo->GetNE();
2169  int meshCounter, spaceCounter, dim = Dimension();
2170 
2171  Array<int> edges;
2172  Array<int> orient;
2173 
2174  v_meshOffsets.SetSize(nv);
2175  e_meshOffsets.SetSize(ne);
2176  f_meshOffsets.SetSize(nf);
2177  p_meshOffsets.SetSize(np);
2178 
2179  v_spaceOffsets.SetSize(nv);
2180  e_spaceOffsets.SetSize(ne);
2181  f_spaceOffsets.SetSize(nf);
2182  p_spaceOffsets.SetSize(np);
2183 
2184  // Get vertex offsets
2185  for (meshCounter = 0; meshCounter < nv; meshCounter++)
2186  {
2187  v_meshOffsets[meshCounter] = meshCounter;
2188  v_spaceOffsets[meshCounter] = meshCounter;
2189  }
2190  spaceCounter = meshCounter;
2191 
2192  // Get edge offsets
2193  for (int e = 0; e < ne; e++)
2194  {
2195  e_meshOffsets[e] = meshCounter;
2196  e_spaceOffsets[e] = spaceCounter;
2197  meshCounter += KnotVec(e)->GetNE() - 1;
2198  spaceCounter += KnotVec(e)->GetNCP() - 2;
2199  }
2200 
2201  // Get face offsets
2202  for (int f = 0; f < nf; f++)
2203  {
2204  f_meshOffsets[f] = meshCounter;
2205  f_spaceOffsets[f] = spaceCounter;
2206 
2207  patchTopo->GetFaceEdges(f, edges, orient);
2208 
2209  meshCounter +=
2210  (KnotVec(edges[0])->GetNE() - 1) *
2211  (KnotVec(edges[1])->GetNE() - 1);
2212  spaceCounter +=
2213  (KnotVec(edges[0])->GetNCP() - 2) *
2214  (KnotVec(edges[1])->GetNCP() - 2);
2215  }
2216 
2217  // Get patch offsets
2218  for (int p = 0; p < np; p++)
2219  {
2220  p_meshOffsets[p] = meshCounter;
2221  p_spaceOffsets[p] = spaceCounter;
2222 
2223  patchTopo->GetElementEdges(p, edges, orient);
2224 
2225  if (dim == 2)
2226  {
2227  meshCounter +=
2228  (KnotVec(edges[0])->GetNE() - 1) *
2229  (KnotVec(edges[1])->GetNE() - 1);
2230  spaceCounter +=
2231  (KnotVec(edges[0])->GetNCP() - 2) *
2232  (KnotVec(edges[1])->GetNCP() - 2);
2233  }
2234  else
2235  {
2236  meshCounter +=
2237  (KnotVec(edges[0])->GetNE() - 1) *
2238  (KnotVec(edges[3])->GetNE() - 1) *
2239  (KnotVec(edges[8])->GetNE() - 1);
2240  spaceCounter +=
2241  (KnotVec(edges[0])->GetNCP() - 2) *
2242  (KnotVec(edges[3])->GetNCP() - 2) *
2243  (KnotVec(edges[8])->GetNCP() - 2);
2244  }
2245  }
2246  NumOfVertices = meshCounter;
2247  NumOfDofs = spaceCounter;
2248 }
2249 
2251 {
2252  int dim = Dimension();
2253  Array<const KnotVector *> kv(dim);
2254 
2255  NumOfElements = 0;
2256  for (int p = 0; p < GetNP(); p++)
2257  {
2258  GetPatchKnotVectors(p, kv);
2259 
2260  int ne = kv[0]->GetNE();
2261  for (int d = 1; d < dim; d++)
2262  {
2263  ne *= kv[d]->GetNE();
2264  }
2265 
2266  NumOfElements += ne;
2267  }
2268 }
2269 
2271 {
2272  int dim = Dimension() - 1;
2273  Array<KnotVector *> kv(dim);
2274 
2275  NumOfBdrElements = 0;
2276  for (int p = 0; p < GetNBP(); p++)
2277  {
2278  GetBdrPatchKnotVectors(p, kv);
2279 
2280  int ne = kv[0]->GetNE();
2281  for (int d = 1; d < dim; d++)
2282  {
2283  ne *= kv[d]->GetNE();
2284  }
2285 
2286  NumOfBdrElements += ne;
2287  }
2288 }
2289 
2291 {
2292  elements.SetSize(GetNE());
2293 
2294  if (Dimension() == 2)
2295  {
2296  Get2DElementTopo(elements);
2297  }
2298  else
2299  {
2300  Get3DElementTopo(elements);
2301  }
2302 }
2303 
2305 {
2306  int el = 0;
2307  int eg = 0;
2308  int ind[4];
2309  NURBSPatchMap p2g(this);
2310  const KnotVector *kv[2];
2311 
2312  for (int p = 0; p < GetNP(); p++)
2313  {
2314  p2g.SetPatchVertexMap(p, kv);
2315  int nx = p2g.nx();
2316  int ny = p2g.ny();
2317 
2318  int patch_attr = patchTopo->GetAttribute(p);
2319 
2320  for (int j = 0; j < ny; j++)
2321  {
2322  for (int i = 0; i < nx; i++)
2323  {
2324  if (activeElem[eg])
2325  {
2326  ind[0] = activeVert[p2g(i, j )];
2327  ind[1] = activeVert[p2g(i+1,j )];
2328  ind[2] = activeVert[p2g(i+1,j+1)];
2329  ind[3] = activeVert[p2g(i, j+1)];
2330 
2331  elements[el] = new Quadrilateral(ind, patch_attr);
2332  el++;
2333  }
2334  eg++;
2335  }
2336  }
2337  }
2338 }
2339 
2341 {
2342  int el = 0;
2343  int eg = 0;
2344  int ind[8];
2345  NURBSPatchMap p2g(this);
2346  const KnotVector *kv[3];
2347 
2348  for (int p = 0; p < GetNP(); p++)
2349  {
2350  p2g.SetPatchVertexMap(p, kv);
2351  int nx = p2g.nx();
2352  int ny = p2g.ny();
2353  int nz = p2g.nz();
2354 
2355  int patch_attr = patchTopo->GetAttribute(p);
2356 
2357  for (int k = 0; k < nz; k++)
2358  {
2359  for (int j = 0; j < ny; j++)
2360  {
2361  for (int i = 0; i < nx; i++)
2362  {
2363  if (activeElem[eg])
2364  {
2365  ind[0] = activeVert[p2g(i, j, k)];
2366  ind[1] = activeVert[p2g(i+1,j, k)];
2367  ind[2] = activeVert[p2g(i+1,j+1,k)];
2368  ind[3] = activeVert[p2g(i, j+1,k)];
2369 
2370  ind[4] = activeVert[p2g(i, j, k+1)];
2371  ind[5] = activeVert[p2g(i+1,j, k+1)];
2372  ind[6] = activeVert[p2g(i+1,j+1,k+1)];
2373  ind[7] = activeVert[p2g(i, j+1,k+1)];
2374 
2375  elements[el] = new Hexahedron(ind, patch_attr);
2376  el++;
2377  }
2378  eg++;
2379  }
2380  }
2381  }
2382  }
2383 }
2384 
2386 {
2387  boundary.SetSize(GetNBE());
2388 
2389  if (Dimension() == 2)
2390  {
2391  Get2DBdrElementTopo(boundary);
2392  }
2393  else
2394  {
2395  Get3DBdrElementTopo(boundary);
2396  }
2397 }
2398 
2400 {
2401  int g_be, l_be;
2402  int ind[2], okv[1];
2403  NURBSPatchMap p2g(this);
2404  const KnotVector *kv[1];
2405 
2406  g_be = l_be = 0;
2407  for (int b = 0; b < GetNBP(); b++)
2408  {
2409  p2g.SetBdrPatchVertexMap(b, kv, okv);
2410  int nx = p2g.nx();
2411 
2412  int bdr_patch_attr = patchTopo->GetBdrAttribute(b);
2413 
2414  for (int i = 0; i < nx; i++)
2415  {
2416  if (activeBdrElem[g_be])
2417  {
2418  int _i = (okv[0] >= 0) ? i : (nx - 1 - i);
2419  ind[0] = activeVert[p2g[_i ]];
2420  ind[1] = activeVert[p2g[_i+1]];
2421 
2422  boundary[l_be] = new Segment(ind, bdr_patch_attr);
2423  l_be++;
2424  }
2425  g_be++;
2426  }
2427  }
2428 }
2429 
2431 {
2432  int g_be, l_be;
2433  int ind[4], okv[2];
2434  NURBSPatchMap p2g(this);
2435  const KnotVector *kv[2];
2436 
2437  g_be = l_be = 0;
2438  for (int b = 0; b < GetNBP(); b++)
2439  {
2440  p2g.SetBdrPatchVertexMap(b, kv, okv);
2441  int nx = p2g.nx();
2442  int ny = p2g.ny();
2443 
2444  int bdr_patch_attr = patchTopo->GetBdrAttribute(b);
2445 
2446  for (int j = 0; j < ny; j++)
2447  {
2448  int _j = (okv[1] >= 0) ? j : (ny - 1 - j);
2449  for (int i = 0; i < nx; i++)
2450  {
2451  if (activeBdrElem[g_be])
2452  {
2453  int _i = (okv[0] >= 0) ? i : (nx - 1 - i);
2454  ind[0] = activeVert[p2g(_i, _j )];
2455  ind[1] = activeVert[p2g(_i+1,_j )];
2456  ind[2] = activeVert[p2g(_i+1,_j+1)];
2457  ind[3] = activeVert[p2g(_i, _j+1)];
2458 
2459  boundary[l_be] = new Quadrilateral(ind, bdr_patch_attr);
2460  l_be++;
2461  }
2462  g_be++;
2463  }
2464  }
2465  }
2466 }
2467 
2469 {
2470  activeDof.SetSize(GetNTotalDof());
2471  activeDof = 0;
2472 
2473  if (Dimension() == 2)
2474  {
2475  Generate2DElementDofTable();
2476  }
2477  else
2478  {
2479  Generate3DElementDofTable();
2480  }
2481 
2482  NumOfActiveDofs = 0;
2483  for (int d = 0; d < GetNTotalDof(); d++)
2484  if (activeDof[d])
2485  {
2486  NumOfActiveDofs++;
2487  activeDof[d] = NumOfActiveDofs;
2488  }
2489 
2490  int *dof = el_dof->GetJ();
2491  int ndof = el_dof->Size_of_connections();
2492  for (int i = 0; i < ndof; i++)
2493  {
2494  dof[i] = activeDof[dof[i]] - 1;
2495  }
2496 }
2497 
2499 {
2500  int el = 0;
2501  int eg = 0;
2502  const KnotVector *kv[2];
2503  NURBSPatchMap p2g(this);
2504 
2505  Array<Connection> el_dof_list;
2506  el_to_patch.SetSize(NumOfActiveElems);
2507  el_to_IJK.SetSize(NumOfActiveElems, 2);
2508 
2509  for (int p = 0; p < GetNP(); p++)
2510  {
2511  p2g.SetPatchDofMap(p, kv);
2512 
2513  // Load dofs
2514  const int ord0 = kv[0]->GetOrder();
2515  const int ord1 = kv[1]->GetOrder();
2516  for (int j = 0; j < kv[1]->GetNKS(); j++)
2517  {
2518  if (kv[1]->isElement(j))
2519  {
2520  for (int i = 0; i < kv[0]->GetNKS(); i++)
2521  {
2522  if (kv[0]->isElement(i))
2523  {
2524  if (activeElem[eg])
2525  {
2526  Connection conn(el,0);
2527  for (int jj = 0; jj <= ord1; jj++)
2528  {
2529  for (int ii = 0; ii <= ord0; ii++)
2530  {
2531  conn.to = DofMap(p2g(i+ii,j+jj));
2532  activeDof[conn.to] = 1;
2533  el_dof_list.Append(conn);
2534  }
2535  }
2536  el_to_patch[el] = p;
2537  el_to_IJK(el,0) = i;
2538  el_to_IJK(el,1) = j;
2539 
2540  el++;
2541  }
2542  eg++;
2543  }
2544  }
2545  }
2546  }
2547  }
2548  // We must NOT sort el_dof_list in this case.
2549  el_dof = new Table(NumOfActiveElems, el_dof_list);
2550 }
2551 
2553 {
2554  int el = 0;
2555  int eg = 0;
2556  const KnotVector *kv[3];
2557  NURBSPatchMap p2g(this);
2558 
2559  Array<Connection> el_dof_list;
2560  el_to_patch.SetSize(NumOfActiveElems);
2561  el_to_IJK.SetSize(NumOfActiveElems, 3);
2562 
2563  for (int p = 0; p < GetNP(); p++)
2564  {
2565  p2g.SetPatchDofMap(p, kv);
2566 
2567  // Load dofs
2568  const int ord0 = kv[0]->GetOrder();
2569  const int ord1 = kv[1]->GetOrder();
2570  const int ord2 = kv[2]->GetOrder();
2571  for (int k = 0; k < kv[2]->GetNKS(); k++)
2572  {
2573  if (kv[2]->isElement(k))
2574  {
2575  for (int j = 0; j < kv[1]->GetNKS(); j++)
2576  {
2577  if (kv[1]->isElement(j))
2578  {
2579  for (int i = 0; i < kv[0]->GetNKS(); i++)
2580  {
2581  if (kv[0]->isElement(i))
2582  {
2583  if (activeElem[eg])
2584  {
2585  Connection conn(el,0);
2586  for (int kk = 0; kk <= ord2; kk++)
2587  {
2588  for (int jj = 0; jj <= ord1; jj++)
2589  {
2590  for (int ii = 0; ii <= ord0; ii++)
2591  {
2592  conn.to = DofMap(p2g(i+ii, j+jj, k+kk));
2593  activeDof[conn.to] = 1;
2594  el_dof_list.Append(conn);
2595  }
2596  }
2597  }
2598 
2599  el_to_patch[el] = p;
2600  el_to_IJK(el,0) = i;
2601  el_to_IJK(el,1) = j;
2602  el_to_IJK(el,2) = k;
2603 
2604  el++;
2605  }
2606  eg++;
2607  }
2608  }
2609  }
2610  }
2611  }
2612  }
2613  }
2614  // We must NOT sort el_dof_list in this case.
2615  el_dof = new Table(NumOfActiveElems, el_dof_list);
2616 }
2617 
2619 {
2620  if (Dimension() == 2)
2621  {
2622  Generate2DBdrElementDofTable();
2623  }
2624  else
2625  {
2626  Generate3DBdrElementDofTable();
2627  }
2628 
2629  int *dof = bel_dof->GetJ();
2630  int ndof = bel_dof->Size_of_connections();
2631  for (int i = 0; i < ndof; i++)
2632  {
2633  dof[i] = activeDof[dof[i]] - 1;
2634  }
2635 }
2636 
2638 {
2639  int gbe = 0;
2640  int lbe = 0, okv[1];
2641  const KnotVector *kv[1];
2642  NURBSPatchMap p2g(this);
2643 
2644  Array<Connection> bel_dof_list;
2645  bel_to_patch.SetSize(NumOfActiveBdrElems);
2646  bel_to_IJK.SetSize(NumOfActiveBdrElems, 1);
2647 
2648  for (int b = 0; b < GetNBP(); b++)
2649  {
2650  p2g.SetBdrPatchDofMap(b, kv, okv);
2651  const int nx = p2g.nx(); // NCP-1
2652 
2653  // Load dofs
2654  const int nks0 = kv[0]->GetNKS();
2655  const int ord0 = kv[0]->GetOrder();
2656  for (int i = 0; i < nks0; i++)
2657  {
2658  if (kv[0]->isElement(i))
2659  {
2660  if (activeBdrElem[gbe])
2661  {
2662  Connection conn(lbe,0);
2663  for (int ii = 0; ii <= ord0; ii++)
2664  {
2665  conn.to = DofMap(p2g[(okv[0] >= 0) ? (i+ii) : (nx-i-ii)]);
2666  bel_dof_list.Append(conn);
2667  }
2668  bel_to_patch[lbe] = b;
2669  bel_to_IJK(lbe,0) = (okv[0] >= 0) ? i : (-1-i);
2670  lbe++;
2671  }
2672  gbe++;
2673  }
2674  }
2675  }
2676  // We must NOT sort bel_dof_list in this case.
2677  bel_dof = new Table(NumOfActiveBdrElems, bel_dof_list);
2678 }
2679 
2681 {
2682  int gbe = 0;
2683  int lbe = 0, okv[2];
2684  const KnotVector *kv[2];
2685  NURBSPatchMap p2g(this);
2686 
2687  Array<Connection> bel_dof_list;
2688  bel_to_patch.SetSize(NumOfActiveBdrElems);
2689  bel_to_IJK.SetSize(NumOfActiveBdrElems, 2);
2690 
2691  for (int b = 0; b < GetNBP(); b++)
2692  {
2693  p2g.SetBdrPatchDofMap(b, kv, okv);
2694  const int nx = p2g.nx(); // NCP0-1
2695  const int ny = p2g.ny(); // NCP1-1
2696 
2697  // Load dofs
2698  const int nks0 = kv[0]->GetNKS();
2699  const int ord0 = kv[0]->GetOrder();
2700  const int nks1 = kv[1]->GetNKS();
2701  const int ord1 = kv[1]->GetOrder();
2702  for (int j = 0; j < nks1; j++)
2703  {
2704  if (kv[1]->isElement(j))
2705  {
2706  for (int i = 0; i < nks0; i++)
2707  {
2708  if (kv[0]->isElement(i))
2709  {
2710  if (activeBdrElem[gbe])
2711  {
2712  Connection conn(lbe,0);
2713  for (int jj = 0; jj <= ord1; jj++)
2714  {
2715  const int _jj = (okv[1] >= 0) ? (j+jj) : (ny-j-jj);
2716  for (int ii = 0; ii <= ord0; ii++)
2717  {
2718  const int _ii = (okv[0] >= 0) ? (i+ii) : (nx-i-ii);
2719  conn.to = DofMap(p2g(_ii, _jj));
2720  bel_dof_list.Append(conn);
2721  }
2722  }
2723  bel_to_patch[lbe] = b;
2724  bel_to_IJK(lbe,0) = (okv[0] >= 0) ? i : (-1-i);
2725  bel_to_IJK(lbe,1) = (okv[1] >= 0) ? j : (-1-j);
2726  lbe++;
2727  }
2728  gbe++;
2729  }
2730  }
2731  }
2732  }
2733  }
2734  // We must NOT sort bel_dof_list in this case.
2735  bel_dof = new Table(NumOfActiveBdrElems, bel_dof_list);
2736 }
2737 
2739 {
2740  lvert_vert.SetSize(GetNV());
2741  for (int gv = 0; gv < GetGNV(); gv++)
2742  if (activeVert[gv] >= 0)
2743  {
2744  lvert_vert[activeVert[gv]] = gv;
2745  }
2746 }
2747 
2749 {
2750  lelem_elem.SetSize(GetNE());
2751  for (int le = 0, ge = 0; ge < GetGNE(); ge++)
2752  if (activeElem[ge])
2753  {
2754  lelem_elem[le++] = ge;
2755  }
2756 }
2757 
2758 void NURBSExtension::LoadFE(int i, const FiniteElement *FE) const
2759 {
2760  const NURBSFiniteElement *NURBSFE =
2761  dynamic_cast<const NURBSFiniteElement *>(FE);
2762 
2763  if (NURBSFE->GetElement() != i)
2764  {
2765  Array<int> dofs;
2766  NURBSFE->SetIJK(el_to_IJK.GetRow(i));
2767  if (el_to_patch[i] != NURBSFE->GetPatch())
2768  {
2769  GetPatchKnotVectors(el_to_patch[i], NURBSFE->KnotVectors());
2770  NURBSFE->SetPatch(el_to_patch[i]);
2771  NURBSFE->SetOrder();
2772  }
2773  el_dof->GetRow(i, dofs);
2774  weights.GetSubVector(dofs, NURBSFE->Weights());
2775  NURBSFE->SetElement(i);
2776  }
2777 }
2778 
2779 void NURBSExtension::LoadBE(int i, const FiniteElement *BE) const
2780 {
2781  const NURBSFiniteElement *NURBSFE =
2782  dynamic_cast<const NURBSFiniteElement *>(BE);
2783 
2784  if (NURBSFE->GetElement() != i)
2785  {
2786  Array<int> dofs;
2787  NURBSFE->SetIJK(bel_to_IJK.GetRow(i));
2788  if (bel_to_patch[i] != NURBSFE->GetPatch())
2789  {
2790  GetBdrPatchKnotVectors(bel_to_patch[i], NURBSFE->KnotVectors());
2791  NURBSFE->SetPatch(bel_to_patch[i]);
2792  NURBSFE->SetOrder();
2793  }
2794  bel_dof->GetRow(i, dofs);
2795  weights.GetSubVector(dofs, NURBSFE->Weights());
2796  NURBSFE->SetElement(i);
2797  }
2798 }
2799 
2801 {
2802  delete el_dof;
2803  delete bel_dof;
2804 
2805  if (patches.Size() == 0)
2806  {
2807  GetPatchNets(Nodes, Dimension());
2808  }
2809 }
2810 
2812 {
2813  if (patches.Size() == 0) { return; }
2814 
2815  SetSolutionVector(Nodes, Dimension());
2816  patches.SetSize(0);
2817 }
2818 
2820 {
2821  if (patches.Size() == 0)
2822  {
2823  mfem_error("NURBSExtension::SetKnotsFromPatches :"
2824  " No patches available!");
2825  }
2826 
2828 
2829  for (int p = 0; p < patches.Size(); p++)
2830  {
2831  GetPatchKnotVectors(p, kv);
2832 
2833  for (int i = 0; i < kv.Size(); i++)
2834  {
2835  *kv[i] = *patches[p]->GetKV(i);
2836  }
2837  }
2838 
2839  SetOrdersFromKnotVectors();
2840 
2841  GenerateOffsets();
2842  CountElements();
2843  CountBdrElements();
2844 
2845  // all elements must be active
2846  NumOfActiveElems = NumOfElements;
2847  activeElem.SetSize(NumOfElements);
2848  activeElem = true;
2849 
2850  GenerateActiveVertices();
2851  InitDofMap();
2852  GenerateElementDofTable();
2853  GenerateActiveBdrElems();
2854  GenerateBdrElementDofTable();
2855 
2856  ConnectBoundaries();
2857 }
2858 
2859 void NURBSExtension::LoadSolution(std::istream &input, GridFunction &sol) const
2860 {
2861  const FiniteElementSpace *fes = sol.FESpace();
2862  MFEM_VERIFY(fes->GetNURBSext() == this, "");
2863 
2864  sol.SetSize(fes->GetVSize());
2865 
2866  Array<const KnotVector *> kv(Dimension());
2867  NURBSPatchMap p2g(this);
2868  const int vdim = fes->GetVDim();
2869 
2870  for (int p = 0; p < GetNP(); p++)
2871  {
2872  skip_comment_lines(input, '#');
2873 
2874  p2g.SetPatchDofMap(p, kv);
2875  const int nx = kv[0]->GetNCP();
2876  const int ny = kv[1]->GetNCP();
2877  const int nz = (kv.Size() == 2) ? 1 : kv[2]->GetNCP();
2878  for (int k = 0; k < nz; k++)
2879  {
2880  for (int j = 0; j < ny; j++)
2881  {
2882  for (int i = 0; i < nx; i++)
2883  {
2884  const int ll = (kv.Size() == 2) ? p2g(i,j) : p2g(i,j,k);
2885  const int l = DofMap(ll);
2886  for (int vd = 0; vd < vdim; vd++)
2887  {
2888  input >> sol(fes->DofToVDof(l,vd));
2889  }
2890  }
2891  }
2892  }
2893  }
2894 }
2895 
2896 void NURBSExtension::PrintSolution(const GridFunction &sol, std::ostream &out)
2897 const
2898 {
2899  const FiniteElementSpace *fes = sol.FESpace();
2900  MFEM_VERIFY(fes->GetNURBSext() == this, "");
2901 
2902  Array<const KnotVector *> kv(Dimension());
2903  NURBSPatchMap p2g(this);
2904  const int vdim = fes->GetVDim();
2905 
2906  for (int p = 0; p < GetNP(); p++)
2907  {
2908  out << "\n# patch " << p << "\n\n";
2909 
2910  p2g.SetPatchDofMap(p, kv);
2911  const int nx = kv[0]->GetNCP();
2912  const int ny = kv[1]->GetNCP();
2913  const int nz = (kv.Size() == 2) ? 1 : kv[2]->GetNCP();
2914  for (int k = 0; k < nz; k++)
2915  {
2916  for (int j = 0; j < ny; j++)
2917  {
2918  for (int i = 0; i < nx; i++)
2919  {
2920  const int ll = (kv.Size() == 2) ? p2g(i,j) : p2g(i,j,k);
2921  const int l = DofMap(ll);
2922  out << sol(fes->DofToVDof(l,0));
2923  for (int vd = 1; vd < vdim; vd++)
2924  {
2925  out << ' ' << sol(fes->DofToVDof(l,vd));
2926  }
2927  out << '\n';
2928  }
2929  }
2930  }
2931  }
2932 }
2933 
2934 void NURBSExtension::DegreeElevate(int rel_degree, int degree)
2935 {
2936  for (int p = 0; p < patches.Size(); p++)
2937  {
2938  for (int dir = 0; dir < patches[p]->GetNKV(); dir++)
2939  {
2940  int oldd = patches[p]->GetKV(dir)->GetOrder();
2941  int newd = std::min(oldd + rel_degree, degree);
2942  if (newd > oldd)
2943  {
2944  patches[p]->DegreeElevate(dir, newd - oldd);
2945  }
2946  }
2947  }
2948 }
2949 
2951 {
2952  for (int p = 0; p < patches.Size(); p++)
2953  {
2954  patches[p]->UniformRefinement();
2955  }
2956 }
2957 
2959 {
2960  Array<int> edges;
2961  Array<int> orient;
2962 
2963  Array<KnotVector *> pkv(Dimension());
2964 
2965  for (int p = 0; p < patches.Size(); p++)
2966  {
2967  patchTopo->GetElementEdges(p, edges, orient);
2968 
2969  if (Dimension()==2)
2970  {
2971  pkv[0] = kv[KnotInd(edges[0])];
2972  pkv[1] = kv[KnotInd(edges[1])];
2973  }
2974  else
2975  {
2976  pkv[0] = kv[KnotInd(edges[0])];
2977  pkv[1] = kv[KnotInd(edges[3])];
2978  pkv[2] = kv[KnotInd(edges[8])];
2979  }
2980 
2981  patches[p]->KnotInsert(pkv);
2982  }
2983 }
2984 
2985 void NURBSExtension::GetPatchNets(const Vector &coords, int vdim)
2986 {
2987  if (Dimension() == 2)
2988  {
2989  Get2DPatchNets(coords, vdim);
2990  }
2991  else
2992  {
2993  Get3DPatchNets(coords, vdim);
2994  }
2995 }
2996 
2997 void NURBSExtension::Get2DPatchNets(const Vector &coords, int vdim)
2998 {
3000  NURBSPatchMap p2g(this);
3001 
3002  patches.SetSize(GetNP());
3003  for (int p = 0; p < GetNP(); p++)
3004  {
3005  p2g.SetPatchDofMap(p, kv);
3006  patches[p] = new NURBSPatch(kv, vdim+1);
3007  NURBSPatch &Patch = *patches[p];
3008 
3009  for (int j = 0; j < kv[1]->GetNCP(); j++)
3010  {
3011  for (int i = 0; i < kv[0]->GetNCP(); i++)
3012  {
3013  const int l = DofMap(p2g(i,j));
3014  for (int d = 0; d < vdim; d++)
3015  {
3016  Patch(i,j,d) = coords(l*vdim + d)*weights(l);
3017  }
3018  Patch(i,j,vdim) = weights(l);
3019  }
3020  }
3021  }
3022 }
3023 
3024 void NURBSExtension::Get3DPatchNets(const Vector &coords, int vdim)
3025 {
3027  NURBSPatchMap p2g(this);
3028 
3029  patches.SetSize(GetNP());
3030  for (int p = 0; p < GetNP(); p++)
3031  {
3032  p2g.SetPatchDofMap(p, kv);
3033  patches[p] = new NURBSPatch(kv, vdim+1);
3034  NURBSPatch &Patch = *patches[p];
3035 
3036  for (int k = 0; k < kv[2]->GetNCP(); k++)
3037  {
3038  for (int j = 0; j < kv[1]->GetNCP(); j++)
3039  {
3040  for (int i = 0; i < kv[0]->GetNCP(); i++)
3041  {
3042  const int l = DofMap(p2g(i,j,k));
3043  for (int d = 0; d < vdim; d++)
3044  {
3045  Patch(i,j,k,d) = coords(l*vdim + d)*weights(l);
3046  }
3047  Patch(i,j,k,vdim) = weights(l);
3048  }
3049  }
3050  }
3051  }
3052 }
3053 
3055 {
3056  if (Dimension() == 2)
3057  {
3058  Set2DSolutionVector(coords, vdim);
3059  }
3060  else
3061  {
3062  Set3DSolutionVector(coords, vdim);
3063  }
3064 }
3065 
3067 {
3069  NURBSPatchMap p2g(this);
3070 
3071  weights.SetSize(GetNDof());
3072  for (int p = 0; p < GetNP(); p++)
3073  {
3074  p2g.SetPatchDofMap(p, kv);
3075  NURBSPatch &Patch = *patches[p];
3076  MFEM_ASSERT(vdim+1 == Patch.GetNC(), "");
3077 
3078  for (int j = 0; j < kv[1]->GetNCP(); j++)
3079  {
3080  for (int i = 0; i < kv[0]->GetNCP(); i++)
3081  {
3082  const int l = p2g(i,j);
3083  for (int d = 0; d < vdim; d++)
3084  {
3085  coords(l*vdim + d) = Patch(i,j,d)/Patch(i,j,vdim);
3086  }
3087  weights(l) = Patch(i,j,vdim);
3088  }
3089  }
3090  delete patches[p];
3091  }
3092 }
3093 
3095 {
3097  NURBSPatchMap p2g(this);
3098 
3099  weights.SetSize(GetNDof());
3100  for (int p = 0; p < GetNP(); p++)
3101  {
3102  p2g.SetPatchDofMap(p, kv);
3103  NURBSPatch &Patch = *patches[p];
3104  MFEM_ASSERT(vdim+1 == Patch.GetNC(), "");
3105 
3106  for (int k = 0; k < kv[2]->GetNCP(); k++)
3107  {
3108  for (int j = 0; j < kv[1]->GetNCP(); j++)
3109  {
3110  for (int i = 0; i < kv[0]->GetNCP(); i++)
3111  {
3112  const int l = p2g(i,j,k);
3113  for (int d = 0; d < vdim; d++)
3114  {
3115  coords(l*vdim + d) = Patch(i,j,k,d)/Patch(i,j,k,vdim);
3116  }
3117  weights(l) = Patch(i,j,k,vdim);
3118  }
3119  }
3120  }
3121  delete patches[p];
3122  }
3123 }
3124 
3125 
3126 #ifdef MFEM_USE_MPI
3128  : NURBSExtension(orig),
3129  partitioning(orig.partitioning ? new int[orig.GetGNE()] : NULL),
3130  gtopo(orig.gtopo),
3131  ldof_group(orig.ldof_group)
3132 {
3133  // Copy the partitioning, if not NULL
3134  if (partitioning)
3135  {
3136  std::memcpy(partitioning, orig.partitioning, orig.GetGNE()*sizeof(int));
3137  }
3138 }
3139 
3141  int *part, const Array<bool> &active_bel)
3142  : gtopo(comm)
3143 {
3144  if (parent->NumOfActiveElems < parent->NumOfElements)
3145  {
3146  // SetActive (BuildGroups?) and the way the weights are copied
3147  // do not support this case
3148  mfem_error("ParNURBSExtension::ParNURBSExtension :\n"
3149  " all elements in the parent must be active!");
3150  }
3151 
3152  patchTopo = parent->patchTopo;
3153  // steal ownership of patchTopo from the 'parent' NURBS extension
3154  if (!parent->own_topo)
3155  {
3156  mfem_error("ParNURBSExtension::ParNURBSExtension :\n"
3157  " parent does not own the patch topology!");
3158  }
3159  own_topo = 1;
3160  parent->own_topo = 0;
3161 
3162  parent->edge_to_knot.Copy(edge_to_knot);
3163 
3164  parent->GetOrders().Copy(mOrders);
3165  mOrder = parent->GetOrder();
3166 
3167  NumOfKnotVectors = parent->GetNKV();
3168  knotVectors.SetSize(NumOfKnotVectors);
3169  for (int i = 0; i < NumOfKnotVectors; i++)
3170  {
3171  knotVectors[i] = new KnotVector(*parent->GetKnotVector(i));
3172  }
3173 
3174  GenerateOffsets();
3175  CountElements();
3176  CountBdrElements();
3177 
3178  // copy 'part' to 'partitioning'
3179  partitioning = new int[GetGNE()];
3180  for (int i = 0; i < GetGNE(); i++)
3181  {
3182  partitioning[i] = part[i];
3183  }
3184  SetActive(partitioning, active_bel);
3185 
3188  // GenerateActiveBdrElems(); // done by SetActive for now
3190 
3191  Table *serial_elem_dof = parent->GetElementDofTable();
3192  BuildGroups(partitioning, *serial_elem_dof);
3193 
3194  weights.SetSize(GetNDof());
3195  // copy weights from parent
3196  for (int gel = 0, lel = 0; gel < GetGNE(); gel++)
3197  {
3198  if (activeElem[gel])
3199  {
3200  int ndofs = el_dof->RowSize(lel);
3201  int *ldofs = el_dof->GetRow(lel);
3202  int *gdofs = serial_elem_dof->GetRow(gel);
3203  for (int i = 0; i < ndofs; i++)
3204  {
3205  weights(ldofs[i]) = parent->weights(gdofs[i]);
3206  }
3207  lel++;
3208  }
3209  }
3210 }
3211 
3213  const ParNURBSExtension *par_parent)
3214  : gtopo(par_parent->gtopo.GetComm())
3215 {
3216  // steal all data from parent
3217  mOrder = parent->mOrder;
3218  Swap(mOrders, parent->mOrders);
3219 
3220  patchTopo = parent->patchTopo;
3221  own_topo = parent->own_topo;
3222  parent->own_topo = 0;
3223 
3224  Swap(edge_to_knot, parent->edge_to_knot);
3225 
3227  Swap(knotVectors, parent->knotVectors);
3228 
3229  NumOfVertices = parent->NumOfVertices;
3230  NumOfElements = parent->NumOfElements;
3232  NumOfDofs = parent->NumOfDofs;
3233 
3234  Swap(v_meshOffsets, parent->v_meshOffsets);
3235  Swap(e_meshOffsets, parent->e_meshOffsets);
3236  Swap(f_meshOffsets, parent->f_meshOffsets);
3237  Swap(p_meshOffsets, parent->p_meshOffsets);
3238 
3243 
3244  Swap(d_to_d, parent->d_to_d);
3245  Swap(master, parent->master);
3246  Swap(slave, parent->slave);
3247 
3251  NumOfActiveDofs = parent->NumOfActiveDofs;
3252 
3253  Swap(activeVert, parent->activeVert);
3254  Swap(activeElem, parent->activeElem);
3255  Swap(activeBdrElem, parent->activeBdrElem);
3256  Swap(activeDof, parent->activeDof);
3257 
3258  el_dof = parent->el_dof;
3259  bel_dof = parent->bel_dof;
3260  parent->el_dof = parent->bel_dof = NULL;
3261 
3262  Swap(el_to_patch, parent->el_to_patch);
3263  Swap(bel_to_patch, parent->bel_to_patch);
3264  Swap(el_to_IJK, parent->el_to_IJK);
3265  Swap(bel_to_IJK, parent->bel_to_IJK);
3266 
3267  Swap(weights, parent->weights);
3268  MFEM_VERIFY(!parent->HavePatches(), "");
3269 
3270  delete parent;
3271 
3272  partitioning = NULL;
3273 
3274  MFEM_VERIFY(par_parent->partitioning,
3275  "parent ParNURBSExtension has no partitioning!");
3276 
3277  // Support for the case when 'parent' is not a local NURBSExtension, i.e.
3278  // NumOfActiveElems is not the same as in 'par_parent'. In that case, we
3279  // assume 'parent' is a global NURBSExtension, i.e. all elements are active.
3280  bool extract_weights = false;
3281  if (NumOfActiveElems != par_parent->NumOfActiveElems)
3282  {
3283  MFEM_ASSERT(NumOfActiveElems == NumOfElements, "internal error");
3284 
3285  SetActive(par_parent->partitioning, par_parent->activeBdrElem);
3287  delete el_dof;
3289  el_to_IJK.DeleteAll();
3291  // GenerateActiveBdrElems(); // done by SetActive for now
3292  delete bel_dof;
3296  extract_weights = true;
3297  }
3298 
3299  Table *glob_elem_dof = GetGlobalElementDofTable();
3300  BuildGroups(par_parent->partitioning, *glob_elem_dof);
3301  if (extract_weights)
3302  {
3303  Vector glob_weights;
3304  Swap(weights, glob_weights);
3305  weights.SetSize(GetNDof());
3306  // Copy the local 'weights' from the 'glob_weights'.
3307  // Assumption: the local element ids follow the global ordering.
3308  for (int gel = 0, lel = 0; gel < GetGNE(); gel++)
3309  {
3310  if (activeElem[gel])
3311  {
3312  int ndofs = el_dof->RowSize(lel);
3313  int *ldofs = el_dof->GetRow(lel);
3314  int *gdofs = glob_elem_dof->GetRow(gel);
3315  for (int i = 0; i < ndofs; i++)
3316  {
3317  weights(ldofs[i]) = glob_weights(gdofs[i]);
3318  }
3319  lel++;
3320  }
3321  }
3322  }
3323  delete glob_elem_dof;
3324 }
3325 
3326 Table *ParNURBSExtension::GetGlobalElementDofTable()
3327 {
3328  if (Dimension() == 2)
3329  {
3330  return Get2DGlobalElementDofTable();
3331  }
3332  else
3333  {
3334  return Get3DGlobalElementDofTable();
3335  }
3336 }
3337 
3338 Table *ParNURBSExtension::Get2DGlobalElementDofTable()
3339 {
3340  int el = 0;
3341  const KnotVector *kv[2];
3342  NURBSPatchMap p2g(this);
3343  Array<Connection> gel_dof_list;
3344 
3345  for (int p = 0; p < GetNP(); p++)
3346  {
3347  p2g.SetPatchDofMap(p, kv);
3348 
3349  // Load dofs
3350  const int ord0 = kv[0]->GetOrder();
3351  const int ord1 = kv[1]->GetOrder();
3352  for (int j = 0; j < kv[1]->GetNKS(); j++)
3353  {
3354  if (kv[1]->isElement(j))
3355  {
3356  for (int i = 0; i < kv[0]->GetNKS(); i++)
3357  {
3358  if (kv[0]->isElement(i))
3359  {
3360  Connection conn(el,0);
3361  for (int jj = 0; jj <= ord1; jj++)
3362  {
3363  for (int ii = 0; ii <= ord0; ii++)
3364  {
3365  conn.to = p2g(i+ii,j+jj);
3366  gel_dof_list.Append(conn);
3367  }
3368  }
3369  el++;
3370  }
3371  }
3372  }
3373  }
3374  }
3375  // We must NOT sort gel_dof_list in this case.
3376  return (new Table(GetGNE(), gel_dof_list));
3377 }
3378 
3379 Table *ParNURBSExtension::Get3DGlobalElementDofTable()
3380 {
3381  int el = 0;
3382  const KnotVector *kv[3];
3383  NURBSPatchMap p2g(this);
3384  Array<Connection> gel_dof_list;
3385 
3386  for (int p = 0; p < GetNP(); p++)
3387  {
3388  p2g.SetPatchDofMap(p, kv);
3389 
3390  // Load dofs
3391  const int ord0 = kv[0]->GetOrder();
3392  const int ord1 = kv[1]->GetOrder();
3393  const int ord2 = kv[2]->GetOrder();
3394  for (int k = 0; k < kv[2]->GetNKS(); k++)
3395  {
3396  if (kv[2]->isElement(k))
3397  {
3398  for (int j = 0; j < kv[1]->GetNKS(); j++)
3399  {
3400  if (kv[1]->isElement(j))
3401  {
3402  for (int i = 0; i < kv[0]->GetNKS(); i++)
3403  {
3404  if (kv[0]->isElement(i))
3405  {
3406  Connection conn(el,0);
3407  for (int kk = 0; kk <= ord2; kk++)
3408  {
3409  for (int jj = 0; jj <= ord1; jj++)
3410  {
3411  for (int ii = 0; ii <= ord0; ii++)
3412  {
3413  conn.to = p2g(i+ii,j+jj,k+kk);
3414  gel_dof_list.Append(conn);
3415  }
3416  }
3417  }
3418  el++;
3419  }
3420  }
3421  }
3422  }
3423  }
3424  }
3425  }
3426  // We must NOT sort gel_dof_list in this case.
3427  return (new Table(GetGNE(), gel_dof_list));
3428 }
3429 
3430 void ParNURBSExtension::SetActive(const int *_partitioning,
3431  const Array<bool> &active_bel)
3432 {
3434  activeElem = false;
3435  NumOfActiveElems = 0;
3436  const int MyRank = gtopo.MyRank();
3437  for (int i = 0; i < GetGNE(); i++)
3438  if (_partitioning[i] == MyRank)
3439  {
3440  activeElem[i] = true;
3441  NumOfActiveElems++;
3442  }
3443 
3444  active_bel.Copy(activeBdrElem);
3445  NumOfActiveBdrElems = 0;
3446  for (int i = 0; i < GetGNBE(); i++)
3447  if (activeBdrElem[i])
3448  {
3450  }
3451 }
3452 
3453 void ParNURBSExtension::BuildGroups(const int *_partitioning,
3454  const Table &elem_dof)
3455 {
3456  Table dof_proc;
3457 
3458  ListOfIntegerSets groups;
3459  IntegerSet group;
3460 
3461  Transpose(elem_dof, dof_proc); // dof_proc is dof_elem
3462  // convert elements to processors
3463  for (int i = 0; i < dof_proc.Size_of_connections(); i++)
3464  {
3465  dof_proc.GetJ()[i] = _partitioning[dof_proc.GetJ()[i]];
3466  }
3467 
3468  // the first group is the local one
3469  int MyRank = gtopo.MyRank();
3470  group.Recreate(1, &MyRank);
3471  groups.Insert(group);
3472 
3473  int dof = 0;
3475  for (int d = 0; d < GetNTotalDof(); d++)
3476  if (activeDof[d])
3477  {
3478  group.Recreate(dof_proc.RowSize(d), dof_proc.GetRow(d));
3479  ldof_group[dof] = groups.Insert(group);
3480 
3481  dof++;
3482  }
3483 
3484  gtopo.Create(groups, 1822);
3485 }
3486 #endif // MFEM_USE_MPI
3487 
3488 
3489 void NURBSPatchMap::GetPatchKnotVectors(int p, const KnotVector *kv[])
3490 {
3491  Ext->patchTopo->GetElementVertices(p, verts);
3492  Ext->patchTopo->GetElementEdges(p, edges, oedge);
3493  if (Ext->Dimension() == 2)
3494  {
3495  kv[0] = Ext->KnotVec(edges[0]);
3496  kv[1] = Ext->KnotVec(edges[1]);
3497  }
3498  else
3499  {
3500  Ext->patchTopo->GetElementFaces(p, faces, oface);
3501 
3502  kv[0] = Ext->KnotVec(edges[0]);
3503  kv[1] = Ext->KnotVec(edges[3]);
3504  kv[2] = Ext->KnotVec(edges[8]);
3505  }
3506  opatch = 0;
3507 }
3508 
3509 void NURBSPatchMap::GetBdrPatchKnotVectors(int p, const KnotVector *kv[],
3510  int *okv)
3511 {
3512  Ext->patchTopo->GetBdrElementVertices(p, verts);
3513  Ext->patchTopo->GetBdrElementEdges(p, edges, oedge);
3514  kv[0] = Ext->KnotVec(edges[0], oedge[0], &okv[0]);
3515  if (Ext->Dimension() == 3)
3516  {
3517  faces.SetSize(1);
3518  Ext->patchTopo->GetBdrElementFace(p, &faces[0], &opatch);
3519  kv[1] = Ext->KnotVec(edges[1], oedge[1], &okv[1]);
3520  }
3521  else
3522  {
3523  opatch = oedge[0];
3524  }
3525 }
3526 
3528 {
3529  GetPatchKnotVectors(p, kv);
3530 
3531  I = kv[0]->GetNE() - 1;
3532  J = kv[1]->GetNE() - 1;
3533 
3534  for (int i = 0; i < verts.Size(); i++)
3535  {
3536  verts[i] = Ext->v_meshOffsets[verts[i]];
3537  }
3538 
3539  for (int i = 0; i < edges.Size(); i++)
3540  {
3541  edges[i] = Ext->e_meshOffsets[edges[i]];
3542  }
3543 
3544  if (Ext->Dimension() == 3)
3545  {
3546  K = kv[2]->GetNE() - 1;
3547 
3548  for (int i = 0; i < faces.Size(); i++)
3549  {
3550  faces[i] = Ext->f_meshOffsets[faces[i]];
3551  }
3552  }
3553 
3554  pOffset = Ext->p_meshOffsets[p];
3555 }
3556 
3558 {
3559  GetPatchKnotVectors(p, kv);
3560 
3561  I = kv[0]->GetNCP() - 2;
3562  J = kv[1]->GetNCP() - 2;
3563 
3564  for (int i = 0; i < verts.Size(); i++)
3565  {
3566  verts[i] = Ext->v_spaceOffsets[verts[i]];
3567  }
3568 
3569  for (int i = 0; i < edges.Size(); i++)
3570  {
3571  edges[i] = Ext->e_spaceOffsets[edges[i]];
3572  }
3573 
3574  if (Ext->Dimension() == 3)
3575  {
3576  K = kv[2]->GetNCP() - 2;
3577 
3578  for (int i = 0; i < faces.Size(); i++)
3579  {
3580  faces[i] = Ext->f_spaceOffsets[faces[i]];
3581  }
3582  }
3583 
3584  pOffset = Ext->p_spaceOffsets[p];
3585 }
3586 
3588  int *okv)
3589 {
3590  GetBdrPatchKnotVectors(p, kv, okv);
3591 
3592  I = kv[0]->GetNE() - 1;
3593 
3594  for (int i = 0; i < verts.Size(); i++)
3595  {
3596  verts[i] = Ext->v_meshOffsets[verts[i]];
3597  }
3598 
3599  if (Ext->Dimension() == 2)
3600  {
3601  pOffset = Ext->e_meshOffsets[edges[0]];
3602  }
3603  else
3604  {
3605  J = kv[1]->GetNE() - 1;
3606 
3607  for (int i = 0; i < edges.Size(); i++)
3608  {
3609  edges[i] = Ext->e_meshOffsets[edges[i]];
3610  }
3611 
3612  pOffset = Ext->f_meshOffsets[faces[0]];
3613  }
3614 }
3615 
3616 void NURBSPatchMap::SetBdrPatchDofMap(int p, const KnotVector *kv[], int *okv)
3617 {
3618  GetBdrPatchKnotVectors(p, kv, okv);
3619 
3620  I = kv[0]->GetNCP() - 2;
3621 
3622  for (int i = 0; i < verts.Size(); i++)
3623  {
3624  verts[i] = Ext->v_spaceOffsets[verts[i]];
3625  }
3626 
3627  if (Ext->Dimension() == 2)
3628  {
3629  pOffset = Ext->e_spaceOffsets[edges[0]];
3630  }
3631  else
3632  {
3633  J = kv[1]->GetNCP() - 2;
3634 
3635  for (int i = 0; i < edges.Size(); i++)
3636  {
3637  edges[i] = Ext->e_spaceOffsets[edges[i]];
3638  }
3639 
3640  pOffset = Ext->f_spaceOffsets[faces[0]];
3641  }
3642 }
3643 
3644 }
Array< KnotVector * > knotVectors
Definition: nurbs.hpp:182
Abstract class for Finite Elements.
Definition: fe.hpp:229
Array< int > f_meshOffsets
Definition: nurbs.hpp:195
KnotVector * DegreeElevate(int t) const
Definition: nurbs.cpp:57
void Create(ListOfIntegerSets &groups, int mpitag)
Array2D< int > bel_to_IJK
Definition: nurbs.hpp:209
int Size() const
Logical size of the array.
Definition: array.hpp:118
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:498
void SetCoordsFromPatches(Vector &Nodes)
Definition: nurbs.cpp:2811
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:347
int GetGNE() const
Definition: nurbs.hpp:346
Array< int > activeVert
Definition: nurbs.hpp:174
void ConnectBoundaries2D(int bnd0, int bnd1)
Definition: nurbs.cpp:1669
KnotVector * KnotVec(int edge)
Definition: nurbs.hpp:543
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:27
int DofToVDof(int dof, int vd, int ndofs=-1) const
Definition: fespace.cpp:143
void GetBdrElementEdges(int i, Array< int > &edges, Array< int > &cor) const
Return the indices and the orientations of all edges of bdr element i.
Definition: mesh.cpp:4064
void DegreeElevate(int dir, int t)
Definition: nurbs.cpp:697
int NumOfElements
Definition: nurbs.hpp:37
NURBSPatch * Interpolate(NURBSPatch &p1, NURBSPatch &p2)
Definition: nurbs.cpp:1073
void LoadBE(int i, const FiniteElement *BE) const
Definition: nurbs.cpp:2779
void Generate2DElementDofTable()
Definition: nurbs.cpp:2498
void MergeGridFunctions(GridFunction *gf_array[], int num_pieces, GridFunction &merged)
Definition: nurbs.cpp:1960
Array< int > ldof_group
Definition: nurbs.hpp:411
int GetPatch() const
Definition: fe.hpp:2867
Array< int > master
Definition: nurbs.hpp:189
Array< const KnotVector * > & KnotVectors() const
Definition: fe.hpp:2871
void Generate3DElementDofTable()
Definition: nurbs.cpp:2552
void KnotInsert(int dir, const KnotVector &knot)
Definition: nurbs.cpp:559
void MergeWeights(Mesh *mesh_array[], int num_pieces)
Definition: nurbs.cpp:1935
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:400
void GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: fespace.cpp:172
void Get3DBdrElementTopo(Array< Element * > &boundary) const
Definition: nurbs.cpp:2430
void UniformRefinement()
Definition: nurbs.cpp:2950
void GetElementFaces(int i, Array< int > &, Array< int > &) const
Return the indices and the orientations of all faces of element i.
Definition: mesh.cpp:4258
int GetNE() const
Definition: nurbs.hpp:48
void Copy(Array &copy) const
Create a copy of the current array.
Definition: array.hpp:795
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Definition: vector.cpp:472
void Rotate3D(double normal[], double angle)
Definition: nurbs.cpp:1024
Array< int > e_meshOffsets
Definition: nurbs.hpp:194
int GetNKV() const
Definition: nurbs.hpp:342
int GetNKS() const
Definition: nurbs.hpp:49
void GetElements()
Count the number of elements.
Definition: nurbs.cpp:100
void SetOrdersFromKnotVectors()
Definition: nurbs.cpp:2153
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
Definition: table.cpp:191
void Generate3DBdrElementDofTable()
Definition: nurbs.cpp:2680
int NumOfControlPoints
Definition: nurbs.hpp:37
Array< int > e_spaceOffsets
Definition: nurbs.hpp:200
void GenerateBdrElementDofTable()
Definition: nurbs.cpp:2618
const KnotVector * GetKnotVector(int i) const
Definition: nurbs.hpp:355
void GetElementLocalToGlobal(Array< int > &lelem_elem)
Definition: nurbs.cpp:2748
void init(int dim_)
Definition: nurbs.cpp:277
Data type quadrilateral element.
Array< int > edge_to_knot
Definition: nurbs.hpp:181
void GenerateActiveVertices()
Definition: nurbs.cpp:1846
Array< bool > activeElem
Definition: nurbs.hpp:175
void GetVertexLocalToGlobal(Array< int > &lvert_vert)
Definition: nurbs.cpp:2738
void Set2DSolutionVector(Vector &Nodes, int vdim)
Definition: nurbs.cpp:3066
Vector knot
Definition: nurbs.hpp:36
int GetElement() const
Definition: fe.hpp:2869
int GetNTotalDof() const
Definition: nurbs.hpp:351
void skip_comment_lines(std::istream &is, const char comment_char)
Definition: text.hpp:27
void DeleteAll()
Delete whole array.
Definition: array.hpp:785
int GetNDof() const
Definition: nurbs.hpp:352
int Dimension() const
Definition: nurbs.hpp:332
Array< int > el_to_patch
Definition: nurbs.hpp:206
void SwapDirections(int dir1, int dir2)
Definition: nurbs.cpp:960
friend class NURBSPatchMap
Definition: nurbs.hpp:162
KnotVector()
Create KnotVector.
Definition: nurbs.hpp:41
static const int MaxOrder
Definition: nurbs.hpp:34
Array< int > v_spaceOffsets
Definition: nurbs.hpp:199
void UniformRefinement(Vector &newknots) const
Definition: nurbs.cpp:86
void Generate2DBdrElementDofTable()
Definition: nurbs.cpp:2637
Data type hexahedron element.
Definition: hexahedron.hpp:22
void ConnectBoundaries3D(int bnd0, int bnd1)
Definition: nurbs.cpp:1745
int dim
Definition: ex3.cpp:48
ParNURBSExtension(const ParNURBSExtension &orig)
Definition: nurbs.cpp:3127
void SetIJK(const int *IJK) const
Definition: fe.hpp:2866
Array< int > d_to_d
Definition: nurbs.hpp:188
void Set3DSolutionVector(Vector &Nodes, int vdim)
Definition: nurbs.cpp:3094
void SetPatch(int p) const
Definition: fe.hpp:2868
void GetPatchKnotVectors(int p, Array< KnotVector * > &kv)
Definition: nurbs.cpp:2059
void GenerateActiveBdrElems()
Definition: nurbs.cpp:1913
int Append(const T &el)
Append element to array, resize if necessary.
Definition: array.hpp:690
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:272
void PrintSolution(const GridFunction &sol, std::ostream &out) const
Definition: nurbs.cpp:2896
void Get3DPatchNets(const Vector &Nodes, int vdim)
Definition: nurbs.cpp:3024
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:146
void Get2DPatchNets(const Vector &Nodes, int vdim)
Definition: nurbs.cpp:2997
bool HavePatches() const
Definition: nurbs.hpp:361
void GenerateOffsets()
Definition: nurbs.cpp:2163
Array< int > activeDof
Definition: nurbs.hpp:177
void GetBdrElementTopo(Array< Element * > &boundary) const
Definition: nurbs.cpp:2385
void Transpose(const Table &A, Table &At, int _ncols_A)
Transpose a Table.
Definition: table.cpp:414
void SetKnotsFromPatches()
Definition: nurbs.cpp:2819
void GetElementTopo(Array< Element * > &elements) const
Definition: nurbs.cpp:2290
const Array< int > & GetOrders() const
Read-only access to the orders of all knot vectors.
Definition: nurbs.hpp:337
Array< int > f_spaceOffsets
Definition: nurbs.hpp:201
Array< bool > activeBdrElem
Definition: nurbs.hpp:176
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:336
int GetNC() const
Definition: nurbs.hpp:127
void LoadSolution(std::istream &input, GridFunction &sol) const
Definition: nurbs.cpp:2859
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:430
void SetElement(int e) const
Definition: fe.hpp:2870
int GetOrder() const
If all orders are identical, return that number. Otherwise, return NURBSFECollection::VariableOrder.
Definition: nurbs.hpp:340
void GetBdrElementFace(int i, int *, int *) const
Return the index and the orientation of the face of bdr element i. (3D)
Definition: mesh.cpp:4296
Array< int > slave
Definition: nurbs.hpp:190
void Swap(Array< T > &, Array< T > &)
Definition: array.hpp:575
void LoadFE(int i, const FiniteElement *FE) const
Definition: nurbs.cpp:2758
const NURBSExtension * GetNURBSext() const
Definition: fespace.hpp:274
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:85
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
Definition: globals.hpp:69
Array< int > v_meshOffsets
Definition: nurbs.hpp:193
int GetNP() const
Definition: nurbs.hpp:333
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:618
void SetBdrPatchVertexMap(int p, const KnotVector *kv[], int *okv)
Definition: nurbs.cpp:3587
void GetBdrPatchKnotVectors(int p, Array< KnotVector * > &kv)
Definition: nurbs.cpp:2100
Array< int > p_meshOffsets
Definition: nurbs.hpp:196
KnotVector & operator=(const KnotVector &kv)
Definition: nurbs.cpp:45
void SetPatchVertexMap(int p, const KnotVector *kv[])
Definition: nurbs.cpp:3527
Helper struct for defining a connectivity table, see Table::MakeFromList.
Definition: table.hpp:27
Table * GetElementDofTable()
Definition: nurbs.hpp:363
Array2D< int > el_to_IJK
Definition: nurbs.hpp:208
void SetOrderFromOrders()
Definition: nurbs.cpp:2139
string direction
void CalcShape(Vector &shape, int i, double xi) const
Definition: nurbs.cpp:132
void CountBdrElements()
Definition: nurbs.cpp:2270
int GetNE() const
Definition: nurbs.hpp:347
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:181
void GetPatchNets(const Vector &Nodes, int vdim)
Definition: nurbs.cpp:2985
void GetElementEdges(int i, Array< int > &edges, Array< int > &cor) const
Return the indices and the orientations of all edges of element i.
Definition: mesh.cpp:4042
NURBSPatch * Revolve3D(NURBSPatch &patch, double n[], double ang, int times)
Definition: nurbs.cpp:1120
Array< int > p_spaceOffsets
Definition: nurbs.hpp:202
void SetSolutionVector(Vector &Nodes, int vdim)
Definition: nurbs.cpp:3054
void Get3DElementTopo(Array< Element * > &elements) const
Definition: nurbs.cpp:2340
void SetBdrPatchDofMap(int p, const KnotVector *kv[], int *okv)
Definition: nurbs.cpp:3616
void ConnectBoundaries()
Definition: nurbs.cpp:1640
int GetNCP() const
Definition: nurbs.hpp:50
void GenerateElementDofTable()
Definition: nurbs.cpp:2468
void Difference(const KnotVector &kv, Vector &diff) const
Definition: nurbs.cpp:243
Vector & Weights() const
Definition: fe.hpp:2872
Array< int > mOrders
Definition: nurbs.hpp:166
virtual void SetOrder() const
Update the NURBSFiniteElement according to the currently set knot vectors.
Definition: fe.hpp:2874
Vector data type.
Definition: vector.hpp:48
void Print(std::ostream &out) const
Definition: nurbs.cpp:125
void DeleteAll()
Delete all dynamically allocated memory, reseting all dimentions to zero.
Definition: array.hpp:390
int GetGNBE() const
Definition: nurbs.hpp:348
int RowSize(int i) const
Definition: table.hpp:108
void CalcDShape(Vector &grad, int i, double xi) const
Definition: nurbs.cpp:158
void DegreeElevate(int rel_degree, int degree=16)
Definition: nurbs.cpp:2934
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:64
void CheckBdrPatches()
Definition: nurbs.cpp:2031
void Get2DElementTopo(Array< Element * > &elements) const
Definition: nurbs.cpp:2304
GroupTopology gtopo
Definition: nurbs.hpp:409
Data type line segment element.
Definition: segment.hpp:22
int GetOrder() const
Definition: nurbs.hpp:51
Array< int > bel_to_patch
Definition: nurbs.hpp:207
void SetPatchDofMap(int p, const KnotVector *kv[])
Definition: nurbs.cpp:3557
void GetBdrElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of boundary element i.
Definition: mesh.hpp:807
void GetElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of element i.
Definition: mesh.hpp:803
void ConvertToPatches(const Vector &Nodes)
Definition: nurbs.cpp:2800
void KnotInsert(Array< KnotVector * > &kv)
Definition: nurbs.cpp:2958
void Get2DBdrElementTopo(Array< Element * > &boundary) const
Definition: nurbs.cpp:2399
int SetLoopDirection(int dir)
Definition: nurbs.cpp:482