MFEM  v3.4
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  v_meshOffsets(orig.v_meshOffsets),
1202  e_meshOffsets(orig.e_meshOffsets),
1203  f_meshOffsets(orig.f_meshOffsets),
1204  p_meshOffsets(orig.p_meshOffsets),
1205  v_spaceOffsets(orig.v_spaceOffsets),
1206  e_spaceOffsets(orig.e_spaceOffsets),
1207  f_spaceOffsets(orig.f_spaceOffsets),
1208  p_spaceOffsets(orig.p_spaceOffsets),
1209  el_dof(orig.el_dof ? new Table(*orig.el_dof) : NULL),
1210  bel_dof(orig.bel_dof ? new Table(*orig.bel_dof) : NULL),
1211  el_to_patch(orig.el_to_patch),
1212  bel_to_patch(orig.bel_to_patch),
1213  el_to_IJK(orig.el_to_IJK),
1214  bel_to_IJK(orig.bel_to_IJK),
1215  patches(orig.patches.Size()) // patches are copied in the body
1216 {
1217  // Copy the knot vectors:
1218  for (int i = 0; i < knotVectors.Size(); i++)
1219  {
1220  knotVectors[i] = new KnotVector(*orig.knotVectors[i]);
1221  }
1222 
1223  // Copy the patches:
1224  for (int p = 0; p < patches.Size(); p++)
1225  {
1226  patches[p] = new NURBSPatch(*orig.patches[p]);
1227  }
1228 }
1229 
1230 NURBSExtension::NURBSExtension(std::istream &input)
1231 {
1232  // Read topology
1233  patchTopo = new Mesh;
1234  patchTopo->LoadPatchTopo(input, edge_to_knot);
1235  own_topo = 1;
1236 
1237  CheckPatches();
1238  // CheckBdrPatches();
1239 
1240  skip_comment_lines(input, '#');
1241 
1242  // Read knotvectors or patches
1243  string ident;
1244  input >> ws >> ident; // 'knotvectors' or 'patches'
1245  if (ident == "knotvectors")
1246  {
1247  input >> NumOfKnotVectors;
1248  knotVectors.SetSize(NumOfKnotVectors);
1249  for (int i = 0; i < NumOfKnotVectors; i++)
1250  {
1251  knotVectors[i] = new KnotVector(input);
1252  }
1253  }
1254  else if (ident == "patches")
1255  {
1256  patches.SetSize(GetNP());
1257  for (int p = 0; p < patches.Size(); p++)
1258  {
1259  skip_comment_lines(input, '#');
1260  patches[p] = new NURBSPatch(input);
1261  }
1262 
1263  NumOfKnotVectors = 0;
1264  for (int i = 0; i < patchTopo->GetNEdges(); i++)
1265  if (NumOfKnotVectors < KnotInd(i))
1266  {
1267  NumOfKnotVectors = KnotInd(i);
1268  }
1269  NumOfKnotVectors++;
1270  knotVectors.SetSize(NumOfKnotVectors);
1271  knotVectors = NULL;
1272 
1273  Array<int> edges, oedge;
1274  for (int p = 0; p < patches.Size(); p++)
1275  {
1276  patchTopo->GetElementEdges(p, edges, oedge);
1277  if (Dimension() == 2)
1278  {
1279  if (knotVectors[KnotInd(edges[0])] == NULL)
1280  {
1281  knotVectors[KnotInd(edges[0])] =
1282  new KnotVector(*patches[p]->GetKV(0));
1283  }
1284  if (knotVectors[KnotInd(edges[1])] == NULL)
1285  {
1286  knotVectors[KnotInd(edges[1])] =
1287  new KnotVector(*patches[p]->GetKV(1));
1288  }
1289  }
1290  else
1291  {
1292  if (knotVectors[KnotInd(edges[0])] == NULL)
1293  {
1294  knotVectors[KnotInd(edges[0])] =
1295  new KnotVector(*patches[p]->GetKV(0));
1296  }
1297  if (knotVectors[KnotInd(edges[3])] == NULL)
1298  {
1299  knotVectors[KnotInd(edges[3])] =
1300  new KnotVector(*patches[p]->GetKV(1));
1301  }
1302  if (knotVectors[KnotInd(edges[8])] == NULL)
1303  {
1304  knotVectors[KnotInd(edges[8])] =
1305  new KnotVector(*patches[p]->GetKV(2));
1306  }
1307  }
1308  }
1309  }
1310  else
1311  {
1312  MFEM_ABORT("invalid section: " << ident);
1313  }
1314 
1315  SetOrdersFromKnotVectors();
1316 
1317  GenerateOffsets();
1318  CountElements();
1319  CountBdrElements();
1320  // NumOfVertices, NumOfElements, NumOfBdrElements, NumOfDofs
1321 
1322  skip_comment_lines(input, '#');
1323 
1324  // Check for a list of mesh elements
1325  if (patches.Size() == 0)
1326  {
1327  input >> ws >> ident;
1328  }
1329  if (patches.Size() == 0 && ident == "mesh_elements")
1330  {
1331  input >> NumOfActiveElems;
1332  activeElem.SetSize(GetGNE());
1333  activeElem = false;
1334  int glob_elem;
1335  for (int i = 0; i < NumOfActiveElems; i++)
1336  {
1337  input >> glob_elem;
1338  activeElem[glob_elem] = true;
1339  }
1340 
1341  skip_comment_lines(input, '#');
1342  input >> ws >> ident;
1343  }
1344  else
1345  {
1346  NumOfActiveElems = NumOfElements;
1347  activeElem.SetSize(NumOfElements);
1348  activeElem = true;
1349  }
1350 
1351  GenerateActiveVertices();
1352  GenerateElementDofTable();
1353  GenerateActiveBdrElems();
1354  GenerateBdrElementDofTable();
1355 
1356  if (patches.Size() == 0)
1357  {
1358  // weights
1359  if (ident == "weights")
1360  {
1361  weights.Load(input, GetNDof());
1362  }
1363  else // e.g. ident = "unitweights" or "autoweights"
1364  {
1365  weights.SetSize(GetNDof());
1366  weights = 1.0;
1367  }
1368  }
1369 }
1370 
1371 NURBSExtension::NURBSExtension(NURBSExtension *parent, int newOrder)
1372 {
1373  patchTopo = parent->patchTopo;
1374  own_topo = 0;
1375 
1376  parent->edge_to_knot.Copy(edge_to_knot);
1377 
1378  NumOfKnotVectors = parent->GetNKV();
1379  knotVectors.SetSize(NumOfKnotVectors);
1380  const Array<int> &pOrders = parent->GetOrders();
1381  for (int i = 0; i < NumOfKnotVectors; i++)
1382  {
1383  if (newOrder > pOrders[i])
1384  {
1385  knotVectors[i] =
1386  parent->GetKnotVector(i)->DegreeElevate(newOrder - pOrders[i]);
1387  }
1388  else
1389  {
1390  knotVectors[i] = new KnotVector(*parent->GetKnotVector(i));
1391  }
1392  }
1393 
1394  // copy some data from parent
1395  NumOfElements = parent->NumOfElements;
1396  NumOfBdrElements = parent->NumOfBdrElements;
1397 
1398  SetOrdersFromKnotVectors();
1399 
1400  GenerateOffsets(); // dof offsets will be different from parent
1401 
1402  NumOfActiveVertices = parent->NumOfActiveVertices;
1403  NumOfActiveElems = parent->NumOfActiveElems;
1404  NumOfActiveBdrElems = parent->NumOfActiveBdrElems;
1405  parent->activeVert.Copy(activeVert);
1406  parent->activeElem.Copy(activeElem);
1407  parent->activeBdrElem.Copy(activeBdrElem);
1408 
1409  GenerateElementDofTable();
1410  GenerateBdrElementDofTable();
1411 
1412  weights.SetSize(GetNDof());
1413  weights = 1.0;
1414 }
1415 
1416 NURBSExtension::NURBSExtension(NURBSExtension *parent,
1417  const Array<int> &newOrders)
1418 {
1419  newOrders.Copy(mOrders);
1420  SetOrderFromOrders();
1421 
1422  patchTopo = parent->patchTopo;
1423  own_topo = 0;
1424 
1425  parent->edge_to_knot.Copy(edge_to_knot);
1426 
1427  NumOfKnotVectors = parent->GetNKV();
1428  MFEM_VERIFY(mOrders.Size() == NumOfKnotVectors, "invalid newOrders array");
1429  knotVectors.SetSize(NumOfKnotVectors);
1430  const Array<int> &pOrders = parent->GetOrders();
1431 
1432  for (int i = 0; i < NumOfKnotVectors; i++)
1433  {
1434  if (mOrders[i] > pOrders[i])
1435  {
1436  knotVectors[i] =
1437  parent->GetKnotVector(i)->DegreeElevate(mOrders[i] - pOrders[i]);
1438  }
1439  else
1440  {
1441  knotVectors[i] = new KnotVector(*parent->GetKnotVector(i));
1442  }
1443  }
1444 
1445  // copy some data from parent
1446  NumOfElements = parent->NumOfElements;
1447  NumOfBdrElements = parent->NumOfBdrElements;
1448 
1449  GenerateOffsets(); // dof offsets will be different from parent
1450 
1451  NumOfActiveVertices = parent->NumOfActiveVertices;
1452  NumOfActiveElems = parent->NumOfActiveElems;
1453  NumOfActiveBdrElems = parent->NumOfActiveBdrElems;
1454  parent->activeVert.Copy(activeVert);
1455  parent->activeElem.Copy(activeElem);
1456  parent->activeBdrElem.Copy(activeBdrElem);
1457 
1458  GenerateElementDofTable();
1459  GenerateBdrElementDofTable();
1460 
1461  weights.SetSize(GetNDof());
1462  weights = 1.0;
1463 }
1464 
1465 NURBSExtension::NURBSExtension(Mesh *mesh_array[], int num_pieces)
1466 {
1467  NURBSExtension *parent = mesh_array[0]->NURBSext;
1468 
1469  if (!parent->own_topo)
1470  {
1471  mfem_error("NURBSExtension::NURBSExtension :\n"
1472  " parent does not own the patch topology!");
1473  }
1474  patchTopo = parent->patchTopo;
1475  own_topo = 1;
1476  parent->own_topo = 0;
1477 
1478  parent->edge_to_knot.Copy(edge_to_knot);
1479 
1480  parent->GetOrders().Copy(mOrders);
1481  mOrder = parent->GetOrder();
1482 
1483  NumOfKnotVectors = parent->GetNKV();
1484  knotVectors.SetSize(NumOfKnotVectors);
1485  for (int i = 0; i < NumOfKnotVectors; i++)
1486  {
1487  knotVectors[i] = new KnotVector(*parent->GetKnotVector(i));
1488  }
1489 
1490  GenerateOffsets();
1491  CountElements();
1492  CountBdrElements();
1493 
1494  // assuming the meshes define a partitioning of all the elements
1495  NumOfActiveElems = NumOfElements;
1496  activeElem.SetSize(NumOfElements);
1497  activeElem = true;
1498 
1499  GenerateActiveVertices();
1500  GenerateElementDofTable();
1501  GenerateActiveBdrElems();
1502  GenerateBdrElementDofTable();
1503 
1504  weights.SetSize(GetNDof());
1505  MergeWeights(mesh_array, num_pieces);
1506 }
1507 
1508 NURBSExtension::~NURBSExtension()
1509 {
1510  if (patches.Size() == 0)
1511  {
1512  delete bel_dof;
1513  delete el_dof;
1514  }
1515 
1516  for (int i = 0; i < knotVectors.Size(); i++)
1517  {
1518  delete knotVectors[i];
1519  }
1520 
1521  for (int i = 0; i < patches.Size(); i++)
1522  {
1523  delete patches[i];
1524  }
1525 
1526  if (own_topo)
1527  {
1528  delete patchTopo;
1529  }
1530 }
1531 
1532 void NURBSExtension::Print(std::ostream &out) const
1533 {
1534  patchTopo->PrintTopo(out, edge_to_knot);
1535  if (patches.Size() == 0)
1536  {
1537  out << "\nknotvectors\n" << NumOfKnotVectors << '\n';
1538  for (int i = 0; i < NumOfKnotVectors; i++)
1539  {
1540  knotVectors[i]->Print(out);
1541  }
1542 
1543  if (NumOfActiveElems < NumOfElements)
1544  {
1545  out << "\nmesh_elements\n" << NumOfActiveElems << '\n';
1546  for (int i = 0; i < NumOfElements; i++)
1547  if (activeElem[i])
1548  {
1549  out << i << '\n';
1550  }
1551  }
1552 
1553  out << "\nweights\n";
1554  weights.Print(out, 1);
1555  }
1556  else
1557  {
1558  out << "\npatches\n";
1559  for (int p = 0; p < patches.Size(); p++)
1560  {
1561  out << "\n# patch " << p << "\n\n";
1562  patches[p]->Print(out);
1563  }
1564  }
1565 }
1566 
1567 void NURBSExtension::PrintCharacteristics(std::ostream &out) const
1568 {
1569  out <<
1570  "NURBS Mesh entity sizes:\n"
1571  "Dimension = " << Dimension() << "\n"
1572  "Unique Orders = ";
1573  Array<int> unique_orders(mOrders);
1574  unique_orders.Sort();
1575  unique_orders.Unique();
1576  unique_orders.Print(out, unique_orders.Size());
1577  out <<
1578  "NumOfKnotVectors = " << GetNKV() << "\n"
1579  "NumOfPatches = " << GetNP() << "\n"
1580  "NumOfBdrPatches = " << GetNBP() << "\n"
1581  "NumOfVertices = " << GetGNV() << "\n"
1582  "NumOfElements = " << GetGNE() << "\n"
1583  "NumOfBdrElements = " << GetGNBE() << "\n"
1584  "NumOfDofs = " << GetNTotalDof() << "\n"
1585  "NumOfActiveVertices = " << GetNV() << "\n"
1586  "NumOfActiveElems = " << GetNE() << "\n"
1587  "NumOfActiveBdrElems = " << GetNBE() << "\n"
1588  "NumOfActiveDofs = " << GetNDof() << '\n';
1589  for (int i = 0; i < NumOfKnotVectors; i++)
1590  {
1591  out << ' ' << i + 1 << ") ";
1592  knotVectors[i]->Print(out);
1593  }
1594  out << endl;
1595 }
1596 
1597 void NURBSExtension::GenerateActiveVertices()
1598 {
1599  int vert[8], nv, g_el, nx, ny, nz, dim = Dimension();
1600 
1601  NURBSPatchMap p2g(this);
1602  const KnotVector *kv[3];
1603 
1604  g_el = 0;
1605  activeVert.SetSize(GetGNV());
1606  activeVert = -1;
1607  for (int p = 0; p < GetNP(); p++)
1608  {
1609  p2g.SetPatchVertexMap(p, kv);
1610 
1611  nx = p2g.nx();
1612  ny = p2g.ny();
1613  nz = (dim == 3) ? p2g.nz() : 1;
1614 
1615  for (int k = 0; k < nz; k++)
1616  {
1617  for (int j = 0; j < ny; j++)
1618  {
1619  for (int i = 0; i < nx; i++)
1620  {
1621  if (activeElem[g_el])
1622  {
1623  if (dim == 2)
1624  {
1625  vert[0] = p2g(i, j );
1626  vert[1] = p2g(i+1,j );
1627  vert[2] = p2g(i+1,j+1);
1628  vert[3] = p2g(i, j+1);
1629  nv = 4;
1630  }
1631  else
1632  {
1633  vert[0] = p2g(i, j, k);
1634  vert[1] = p2g(i+1,j, k);
1635  vert[2] = p2g(i+1,j+1,k);
1636  vert[3] = p2g(i, j+1,k);
1637 
1638  vert[4] = p2g(i, j, k+1);
1639  vert[5] = p2g(i+1,j, k+1);
1640  vert[6] = p2g(i+1,j+1,k+1);
1641  vert[7] = p2g(i, j+1,k+1);
1642  nv = 8;
1643  }
1644 
1645  for (int v = 0; v < nv; v++)
1646  {
1647  activeVert[vert[v]] = 1;
1648  }
1649  }
1650  g_el++;
1651  }
1652  }
1653  }
1654  }
1655 
1656  NumOfActiveVertices = 0;
1657  for (int i = 0; i < GetGNV(); i++)
1658  if (activeVert[i] == 1)
1659  {
1660  activeVert[i] = NumOfActiveVertices++;
1661  }
1662 }
1663 
1664 void NURBSExtension::GenerateActiveBdrElems()
1665 {
1666  int dim = Dimension();
1667  Array<KnotVector *> kv(dim);
1668 
1669  activeBdrElem.SetSize(GetGNBE());
1670  if (GetGNE() == GetNE())
1671  {
1672  activeBdrElem = true;
1673  NumOfActiveBdrElems = GetGNBE();
1674  return;
1675  }
1676  activeBdrElem = false;
1677  NumOfActiveBdrElems = 0;
1678  // the mesh will generate the actual boundary including boundary
1679  // elements that are not on boundary patches. we use this for
1680  // visualization of processor boundaries
1681 
1682  // TODO: generate actual boundary?
1683 }
1684 
1685 
1686 void NURBSExtension::MergeWeights(Mesh *mesh_array[], int num_pieces)
1687 {
1688  Array<int> lelem_elem;
1689 
1690  for (int i = 0; i < num_pieces; i++)
1691  {
1692  NURBSExtension *lext = mesh_array[i]->NURBSext;
1693 
1694  lext->GetElementLocalToGlobal(lelem_elem);
1695 
1696  for (int lel = 0; lel < lext->GetNE(); lel++)
1697  {
1698  int gel = lelem_elem[lel];
1699 
1700  int nd = el_dof->RowSize(gel);
1701  int *gdofs = el_dof->GetRow(gel);
1702  int *ldofs = lext->el_dof->GetRow(lel);
1703  for (int j = 0; j < nd; j++)
1704  {
1705  weights(gdofs[j]) = lext->weights(ldofs[j]);
1706  }
1707  }
1708  }
1709 }
1710 
1711 void NURBSExtension::MergeGridFunctions(
1712  GridFunction *gf_array[], int num_pieces, GridFunction &merged)
1713 {
1714  FiniteElementSpace *gfes = merged.FESpace();
1715  Array<int> lelem_elem, dofs;
1716  Vector lvec;
1717 
1718  for (int i = 0; i < num_pieces; i++)
1719  {
1720  FiniteElementSpace *lfes = gf_array[i]->FESpace();
1721  NURBSExtension *lext = lfes->GetMesh()->NURBSext;
1722 
1723  lext->GetElementLocalToGlobal(lelem_elem);
1724 
1725  for (int lel = 0; lel < lext->GetNE(); lel++)
1726  {
1727  lfes->GetElementVDofs(lel, dofs);
1728  gf_array[i]->GetSubVector(dofs, lvec);
1729 
1730  gfes->GetElementVDofs(lelem_elem[lel], dofs);
1731  merged.SetSubVector(dofs, lvec);
1732  }
1733  }
1734 }
1735 
1736 void NURBSExtension::CheckPatches()
1737 {
1738  Array<int> edges;
1739  Array<int> oedge;
1740 
1741  for (int p = 0; p < GetNP(); p++)
1742  {
1743  patchTopo->GetElementEdges(p, edges, oedge);
1744 
1745  for (int i = 0; i < edges.Size(); i++)
1746  {
1747  edges[i] = edge_to_knot[edges[i]];
1748  if (oedge[i] < 0)
1749  {
1750  edges[i] = -1 - edges[i];
1751  }
1752  }
1753 
1754  if ((Dimension() == 2 &&
1755  (edges[0] != -1 - edges[2] || edges[1] != -1 - edges[3])) ||
1756 
1757  (Dimension() == 3 &&
1758  (edges[0] != edges[2] || edges[0] != edges[4] ||
1759  edges[0] != edges[6] || edges[1] != edges[3] ||
1760  edges[1] != edges[5] || edges[1] != edges[7] ||
1761  edges[8] != edges[9] || edges[8] != edges[10] ||
1762  edges[8] != edges[11])))
1763  {
1764  mfem::err << "NURBSExtension::CheckPatch (patch = " << p
1765  << ")\n Inconsistent edge-to-knot mapping!\n";
1766  mfem_error();
1767  }
1768 
1769  if ((Dimension() == 2 &&
1770  (edges[0] < 0 || edges[1] < 0)) ||
1771 
1772  (Dimension() == 3 &&
1773  (edges[0] < 0 || edges[3] < 0 || edges[8] < 0)))
1774  {
1775  mfem::err << "NURBSExtension::CheckPatch (patch = " << p
1776  << ") : Bad orientation!\n";
1777  mfem_error();
1778  }
1779  }
1780 }
1781 
1782 void NURBSExtension::CheckBdrPatches()
1783 {
1784  Array<int> edges;
1785  Array<int> oedge;
1786 
1787  for (int p = 0; p < GetNBP(); p++)
1788  {
1789  patchTopo->GetBdrElementEdges(p, edges, oedge);
1790 
1791  for (int i = 0; i < edges.Size(); i++)
1792  {
1793  edges[i] = edge_to_knot[edges[i]];
1794  if (oedge[i] < 0)
1795  {
1796  edges[i] = -1 - edges[i];
1797  }
1798  }
1799 
1800  if ((Dimension() == 2 && (edges[0] < 0)) ||
1801  (Dimension() == 3 && (edges[0] < 0 || edges[1] < 0)))
1802  {
1803  mfem::err << "NURBSExtension::CheckBdrPatch (boundary patch = "
1804  << p << ") : Bad orientation!\n";
1805  mfem_error();
1806  }
1807  }
1808 }
1809 
1810 void NURBSExtension::GetPatchKnotVectors(int p, Array<KnotVector *> &kv)
1811 {
1812  Array<int> edges, orient;
1813 
1814  kv.SetSize(Dimension());
1815  patchTopo->GetElementEdges(p, edges, orient);
1816 
1817  if (Dimension() == 2)
1818  {
1819  kv[0] = KnotVec(edges[0]);
1820  kv[1] = KnotVec(edges[1]);
1821  }
1822  else
1823  {
1824  kv[0] = KnotVec(edges[0]);
1825  kv[1] = KnotVec(edges[3]);
1826  kv[2] = KnotVec(edges[8]);
1827  }
1828 }
1829 
1830 void NURBSExtension::GetPatchKnotVectors(int p, Array<const KnotVector *> &kv)
1831 const
1832 {
1833  Array<int> edges, orient;
1834 
1835  kv.SetSize(Dimension());
1836  patchTopo->GetElementEdges(p, edges, orient);
1837 
1838  if (Dimension() == 2)
1839  {
1840  kv[0] = KnotVec(edges[0]);
1841  kv[1] = KnotVec(edges[1]);
1842  }
1843  else
1844  {
1845  kv[0] = KnotVec(edges[0]);
1846  kv[1] = KnotVec(edges[3]);
1847  kv[2] = KnotVec(edges[8]);
1848  }
1849 }
1850 
1851 void NURBSExtension::GetBdrPatchKnotVectors(int p, Array<KnotVector *> &kv)
1852 {
1853  Array<int> edges;
1854  Array<int> orient;
1855 
1856  kv.SetSize(Dimension() - 1);
1857  patchTopo->GetBdrElementEdges(p, edges, orient);
1858 
1859  if (Dimension() == 2)
1860  {
1861  kv[0] = KnotVec(edges[0]);
1862  }
1863  else
1864  {
1865  kv[0] = KnotVec(edges[0]);
1866  kv[1] = KnotVec(edges[1]);
1867  }
1868 }
1869 
1870 void NURBSExtension::GetBdrPatchKnotVectors(
1871  int p, Array<const KnotVector *> &kv) const
1872 {
1873  Array<int> edges;
1874  Array<int> orient;
1875 
1876  kv.SetSize(Dimension() - 1);
1877  patchTopo->GetBdrElementEdges(p, edges, orient);
1878 
1879  if (Dimension() == 2)
1880  {
1881  kv[0] = KnotVec(edges[0]);
1882  }
1883  else
1884  {
1885  kv[0] = KnotVec(edges[0]);
1886  kv[1] = KnotVec(edges[1]);
1887  }
1888 }
1889 
1890 void NURBSExtension::SetOrderFromOrders()
1891 {
1892  MFEM_VERIFY(mOrders.Size() > 0, "");
1893  mOrder = mOrders[0];
1894  for (int i = 1; i < mOrders.Size(); i++)
1895  {
1896  if (mOrders[i] != mOrder)
1897  {
1898  mOrder = NURBSFECollection::VariableOrder;
1899  return;
1900  }
1901  }
1902 }
1903 
1904 void NURBSExtension::SetOrdersFromKnotVectors()
1905 {
1906  mOrders.SetSize(NumOfKnotVectors);
1907  for (int i = 0; i < NumOfKnotVectors; i++)
1908  {
1909  mOrders[i] = knotVectors[i]->GetOrder();
1910  }
1911  SetOrderFromOrders();
1912 }
1913 
1914 void NURBSExtension::GenerateOffsets()
1915 {
1916  int nv = patchTopo->GetNV();
1917  int ne = patchTopo->GetNEdges();
1918  int nf = patchTopo->GetNFaces();
1919  int np = patchTopo->GetNE();
1920  int meshCounter, spaceCounter, dim = Dimension();
1921 
1922  Array<int> edges;
1923  Array<int> orient;
1924 
1925  v_meshOffsets.SetSize(nv);
1926  e_meshOffsets.SetSize(ne);
1927  f_meshOffsets.SetSize(nf);
1928  p_meshOffsets.SetSize(np);
1929 
1930  v_spaceOffsets.SetSize(nv);
1931  e_spaceOffsets.SetSize(ne);
1932  f_spaceOffsets.SetSize(nf);
1933  p_spaceOffsets.SetSize(np);
1934 
1935  // Get vertex offsets
1936  for (meshCounter = 0; meshCounter < nv; meshCounter++)
1937  {
1938  v_meshOffsets[meshCounter] = meshCounter;
1939  v_spaceOffsets[meshCounter] = meshCounter;
1940  }
1941  spaceCounter = meshCounter;
1942 
1943  // Get edge offsets
1944  for (int e = 0; e < ne; e++)
1945  {
1946  e_meshOffsets[e] = meshCounter;
1947  e_spaceOffsets[e] = spaceCounter;
1948  meshCounter += KnotVec(e)->GetNE() - 1;
1949  spaceCounter += KnotVec(e)->GetNCP() - 2;
1950  }
1951 
1952  // Get face offsets
1953  for (int f = 0; f < nf; f++)
1954  {
1955  f_meshOffsets[f] = meshCounter;
1956  f_spaceOffsets[f] = spaceCounter;
1957 
1958  patchTopo->GetFaceEdges(f, edges, orient);
1959 
1960  meshCounter +=
1961  (KnotVec(edges[0])->GetNE() - 1) *
1962  (KnotVec(edges[1])->GetNE() - 1);
1963  spaceCounter +=
1964  (KnotVec(edges[0])->GetNCP() - 2) *
1965  (KnotVec(edges[1])->GetNCP() - 2);
1966  }
1967 
1968  // Get patch offsets
1969  for (int p = 0; p < np; p++)
1970  {
1971  p_meshOffsets[p] = meshCounter;
1972  p_spaceOffsets[p] = spaceCounter;
1973 
1974  patchTopo->GetElementEdges(p, edges, orient);
1975 
1976  if (dim == 2)
1977  {
1978  meshCounter +=
1979  (KnotVec(edges[0])->GetNE() - 1) *
1980  (KnotVec(edges[1])->GetNE() - 1);
1981  spaceCounter +=
1982  (KnotVec(edges[0])->GetNCP() - 2) *
1983  (KnotVec(edges[1])->GetNCP() - 2);
1984  }
1985  else
1986  {
1987  meshCounter +=
1988  (KnotVec(edges[0])->GetNE() - 1) *
1989  (KnotVec(edges[3])->GetNE() - 1) *
1990  (KnotVec(edges[8])->GetNE() - 1);
1991  spaceCounter +=
1992  (KnotVec(edges[0])->GetNCP() - 2) *
1993  (KnotVec(edges[3])->GetNCP() - 2) *
1994  (KnotVec(edges[8])->GetNCP() - 2);
1995  }
1996  }
1997  NumOfVertices = meshCounter;
1998  NumOfDofs = spaceCounter;
1999 }
2000 
2001 void NURBSExtension::CountElements()
2002 {
2003  int dim = Dimension();
2004  Array<const KnotVector *> kv(dim);
2005 
2006  NumOfElements = 0;
2007  for (int p = 0; p < GetNP(); p++)
2008  {
2009  GetPatchKnotVectors(p, kv);
2010 
2011  int ne = kv[0]->GetNE();
2012  for (int d = 1; d < dim; d++)
2013  {
2014  ne *= kv[d]->GetNE();
2015  }
2016 
2017  NumOfElements += ne;
2018  }
2019 }
2020 
2021 void NURBSExtension::CountBdrElements()
2022 {
2023  int dim = Dimension() - 1;
2024  Array<KnotVector *> kv(dim);
2025 
2026  NumOfBdrElements = 0;
2027  for (int p = 0; p < GetNBP(); p++)
2028  {
2029  GetBdrPatchKnotVectors(p, kv);
2030 
2031  int ne = kv[0]->GetNE();
2032  for (int d = 1; d < dim; d++)
2033  {
2034  ne *= kv[d]->GetNE();
2035  }
2036 
2037  NumOfBdrElements += ne;
2038  }
2039 }
2040 
2041 void NURBSExtension::GetElementTopo(Array<Element *> &elements) const
2042 {
2043  elements.SetSize(GetNE());
2044 
2045  if (Dimension() == 2)
2046  {
2047  Get2DElementTopo(elements);
2048  }
2049  else
2050  {
2051  Get3DElementTopo(elements);
2052  }
2053 }
2054 
2055 void NURBSExtension::Get2DElementTopo(Array<Element *> &elements) const
2056 {
2057  int el = 0;
2058  int eg = 0;
2059  int ind[4];
2060  NURBSPatchMap p2g(this);
2061  const KnotVector *kv[2];
2062 
2063  for (int p = 0; p < GetNP(); p++)
2064  {
2065  p2g.SetPatchVertexMap(p, kv);
2066  int nx = p2g.nx();
2067  int ny = p2g.ny();
2068 
2069  int patch_attr = patchTopo->GetAttribute(p);
2070 
2071  for (int j = 0; j < ny; j++)
2072  {
2073  for (int i = 0; i < nx; i++)
2074  {
2075  if (activeElem[eg])
2076  {
2077  ind[0] = activeVert[p2g(i, j )];
2078  ind[1] = activeVert[p2g(i+1,j )];
2079  ind[2] = activeVert[p2g(i+1,j+1)];
2080  ind[3] = activeVert[p2g(i, j+1)];
2081 
2082  elements[el] = new Quadrilateral(ind, patch_attr);
2083  el++;
2084  }
2085  eg++;
2086  }
2087  }
2088  }
2089 }
2090 
2091 void NURBSExtension::Get3DElementTopo(Array<Element *> &elements) const
2092 {
2093  int el = 0;
2094  int eg = 0;
2095  int ind[8];
2096  NURBSPatchMap p2g(this);
2097  const KnotVector *kv[3];
2098 
2099  for (int p = 0; p < GetNP(); p++)
2100  {
2101  p2g.SetPatchVertexMap(p, kv);
2102  int nx = p2g.nx();
2103  int ny = p2g.ny();
2104  int nz = p2g.nz();
2105 
2106  int patch_attr = patchTopo->GetAttribute(p);
2107 
2108  for (int k = 0; k < nz; k++)
2109  {
2110  for (int j = 0; j < ny; j++)
2111  {
2112  for (int i = 0; i < nx; i++)
2113  {
2114  if (activeElem[eg])
2115  {
2116  ind[0] = activeVert[p2g(i, j, k)];
2117  ind[1] = activeVert[p2g(i+1,j, k)];
2118  ind[2] = activeVert[p2g(i+1,j+1,k)];
2119  ind[3] = activeVert[p2g(i, j+1,k)];
2120 
2121  ind[4] = activeVert[p2g(i, j, k+1)];
2122  ind[5] = activeVert[p2g(i+1,j, k+1)];
2123  ind[6] = activeVert[p2g(i+1,j+1,k+1)];
2124  ind[7] = activeVert[p2g(i, j+1,k+1)];
2125 
2126  elements[el] = new Hexahedron(ind, patch_attr);
2127  el++;
2128  }
2129  eg++;
2130  }
2131  }
2132  }
2133  }
2134 }
2135 
2136 void NURBSExtension::GetBdrElementTopo(Array<Element *> &boundary) const
2137 {
2138  boundary.SetSize(GetNBE());
2139 
2140  if (Dimension() == 2)
2141  {
2142  Get2DBdrElementTopo(boundary);
2143  }
2144  else
2145  {
2146  Get3DBdrElementTopo(boundary);
2147  }
2148 }
2149 
2150 void NURBSExtension::Get2DBdrElementTopo(Array<Element *> &boundary) const
2151 {
2152  int g_be, l_be;
2153  int ind[2], okv[1];
2154  NURBSPatchMap p2g(this);
2155  const KnotVector *kv[1];
2156 
2157  g_be = l_be = 0;
2158  for (int b = 0; b < GetNBP(); b++)
2159  {
2160  p2g.SetBdrPatchVertexMap(b, kv, okv);
2161  int nx = p2g.nx();
2162 
2163  int bdr_patch_attr = patchTopo->GetBdrAttribute(b);
2164 
2165  for (int i = 0; i < nx; i++)
2166  {
2167  if (activeBdrElem[g_be])
2168  {
2169  int _i = (okv[0] >= 0) ? i : (nx - 1 - i);
2170  ind[0] = activeVert[p2g[_i ]];
2171  ind[1] = activeVert[p2g[_i+1]];
2172 
2173  boundary[l_be] = new Segment(ind, bdr_patch_attr);
2174  l_be++;
2175  }
2176  g_be++;
2177  }
2178  }
2179 }
2180 
2181 void NURBSExtension::Get3DBdrElementTopo(Array<Element *> &boundary) const
2182 {
2183  int g_be, l_be;
2184  int ind[4], okv[2];
2185  NURBSPatchMap p2g(this);
2186  const KnotVector *kv[2];
2187 
2188  g_be = l_be = 0;
2189  for (int b = 0; b < GetNBP(); b++)
2190  {
2191  p2g.SetBdrPatchVertexMap(b, kv, okv);
2192  int nx = p2g.nx();
2193  int ny = p2g.ny();
2194 
2195  int bdr_patch_attr = patchTopo->GetBdrAttribute(b);
2196 
2197  for (int j = 0; j < ny; j++)
2198  {
2199  int _j = (okv[1] >= 0) ? j : (ny - 1 - j);
2200  for (int i = 0; i < nx; i++)
2201  {
2202  if (activeBdrElem[g_be])
2203  {
2204  int _i = (okv[0] >= 0) ? i : (nx - 1 - i);
2205  ind[0] = activeVert[p2g(_i, _j )];
2206  ind[1] = activeVert[p2g(_i+1,_j )];
2207  ind[2] = activeVert[p2g(_i+1,_j+1)];
2208  ind[3] = activeVert[p2g(_i, _j+1)];
2209 
2210  boundary[l_be] = new Quadrilateral(ind, bdr_patch_attr);
2211  l_be++;
2212  }
2213  g_be++;
2214  }
2215  }
2216  }
2217 }
2218 
2219 void NURBSExtension::GenerateElementDofTable()
2220 {
2221  activeDof.SetSize(GetNTotalDof());
2222  activeDof = 0;
2223 
2224  if (Dimension() == 2)
2225  {
2226  Generate2DElementDofTable();
2227  }
2228  else
2229  {
2230  Generate3DElementDofTable();
2231  }
2232 
2233  NumOfActiveDofs = 0;
2234  for (int d = 0; d < GetNTotalDof(); d++)
2235  if (activeDof[d])
2236  {
2237  NumOfActiveDofs++;
2238  activeDof[d] = NumOfActiveDofs;
2239  }
2240 
2241  int *dof = el_dof->GetJ();
2242  int ndof = el_dof->Size_of_connections();
2243  for (int i = 0; i < ndof; i++)
2244  {
2245  dof[i] = activeDof[dof[i]] - 1;
2246  }
2247 }
2248 
2249 void NURBSExtension::Generate2DElementDofTable()
2250 {
2251  int el = 0;
2252  int eg = 0;
2253  const KnotVector *kv[2];
2254  NURBSPatchMap p2g(this);
2255 
2256  Array<Connection> el_dof_list;
2257  el_to_patch.SetSize(NumOfActiveElems);
2258  el_to_IJK.SetSize(NumOfActiveElems, 2);
2259 
2260  for (int p = 0; p < GetNP(); p++)
2261  {
2262  p2g.SetPatchDofMap(p, kv);
2263 
2264  // Load dofs
2265  const int ord0 = kv[0]->GetOrder();
2266  const int ord1 = kv[1]->GetOrder();
2267  for (int j = 0; j < kv[1]->GetNKS(); j++)
2268  {
2269  if (kv[1]->isElement(j))
2270  {
2271  for (int i = 0; i < kv[0]->GetNKS(); i++)
2272  {
2273  if (kv[0]->isElement(i))
2274  {
2275  if (activeElem[eg])
2276  {
2277  Connection conn(el,0);
2278  for (int jj = 0; jj <= ord1; jj++)
2279  {
2280  for (int ii = 0; ii <= ord0; ii++)
2281  {
2282  conn.to = p2g(i+ii,j+jj);
2283  activeDof[conn.to] = 1;
2284  el_dof_list.Append(conn);
2285  }
2286  }
2287  el_to_patch[el] = p;
2288  el_to_IJK(el,0) = i;
2289  el_to_IJK(el,1) = j;
2290 
2291  el++;
2292  }
2293  eg++;
2294  }
2295  }
2296  }
2297  }
2298  }
2299  // We must NOT sort el_dof_list in this case.
2300  el_dof = new Table(NumOfActiveElems, el_dof_list);
2301 }
2302 
2303 void NURBSExtension::Generate3DElementDofTable()
2304 {
2305  int el = 0;
2306  int eg = 0;
2307  const KnotVector *kv[3];
2308  NURBSPatchMap p2g(this);
2309 
2310  Array<Connection> el_dof_list;
2311  el_to_patch.SetSize(NumOfActiveElems);
2312  el_to_IJK.SetSize(NumOfActiveElems, 3);
2313 
2314  for (int p = 0; p < GetNP(); p++)
2315  {
2316  p2g.SetPatchDofMap(p, kv);
2317 
2318  // Load dofs
2319  const int ord0 = kv[0]->GetOrder();
2320  const int ord1 = kv[1]->GetOrder();
2321  const int ord2 = kv[2]->GetOrder();
2322  for (int k = 0; k < kv[2]->GetNKS(); k++)
2323  {
2324  if (kv[2]->isElement(k))
2325  {
2326  for (int j = 0; j < kv[1]->GetNKS(); j++)
2327  {
2328  if (kv[1]->isElement(j))
2329  {
2330  for (int i = 0; i < kv[0]->GetNKS(); i++)
2331  {
2332  if (kv[0]->isElement(i))
2333  {
2334  if (activeElem[eg])
2335  {
2336  Connection conn(el,0);
2337  for (int kk = 0; kk <= ord2; kk++)
2338  {
2339  for (int jj = 0; jj <= ord1; jj++)
2340  {
2341  for (int ii = 0; ii <= ord0; ii++)
2342  {
2343  conn.to = p2g(i+ii, j+jj, k+kk);
2344  activeDof[conn.to] = 1;
2345  el_dof_list.Append(conn);
2346  }
2347  }
2348  }
2349 
2350  el_to_patch[el] = p;
2351  el_to_IJK(el,0) = i;
2352  el_to_IJK(el,1) = j;
2353  el_to_IJK(el,2) = k;
2354 
2355  el++;
2356  }
2357  eg++;
2358  }
2359  }
2360  }
2361  }
2362  }
2363  }
2364  }
2365  // We must NOT sort el_dof_list in this case.
2366  el_dof = new Table(NumOfActiveElems, el_dof_list);
2367 }
2368 
2369 void NURBSExtension::GenerateBdrElementDofTable()
2370 {
2371  if (Dimension() == 2)
2372  {
2373  Generate2DBdrElementDofTable();
2374  }
2375  else
2376  {
2377  Generate3DBdrElementDofTable();
2378  }
2379 
2380  int *dof = bel_dof->GetJ();
2381  int ndof = bel_dof->Size_of_connections();
2382  for (int i = 0; i < ndof; i++)
2383  {
2384  dof[i] = activeDof[dof[i]] - 1;
2385  }
2386 }
2387 
2388 void NURBSExtension::Generate2DBdrElementDofTable()
2389 {
2390  int gbe = 0;
2391  int lbe = 0, okv[1];
2392  const KnotVector *kv[1];
2393  NURBSPatchMap p2g(this);
2394 
2395  Array<Connection> bel_dof_list;
2396  bel_to_patch.SetSize(NumOfActiveBdrElems);
2397  bel_to_IJK.SetSize(NumOfActiveBdrElems, 1);
2398 
2399  for (int b = 0; b < GetNBP(); b++)
2400  {
2401  p2g.SetBdrPatchDofMap(b, kv, okv);
2402  const int nx = p2g.nx(); // NCP-1
2403 
2404  // Load dofs
2405  const int nks0 = kv[0]->GetNKS();
2406  const int ord0 = kv[0]->GetOrder();
2407  for (int i = 0; i < nks0; i++)
2408  {
2409  if (kv[0]->isElement(i))
2410  {
2411  if (activeBdrElem[gbe])
2412  {
2413  Connection conn(lbe,0);
2414  for (int ii = 0; ii <= ord0; ii++)
2415  {
2416  conn.to = p2g[(okv[0] >= 0) ? (i+ii) : (nx-i-ii)];
2417  bel_dof_list.Append(conn);
2418  }
2419  bel_to_patch[lbe] = b;
2420  bel_to_IJK(lbe,0) = (okv[0] >= 0) ? i : (-1-i);
2421  lbe++;
2422  }
2423  gbe++;
2424  }
2425  }
2426  }
2427  // We must NOT sort bel_dof_list in this case.
2428  bel_dof = new Table(NumOfActiveBdrElems, bel_dof_list);
2429 }
2430 
2431 void NURBSExtension::Generate3DBdrElementDofTable()
2432 {
2433  int gbe = 0;
2434  int lbe = 0, okv[2];
2435  const KnotVector *kv[2];
2436  NURBSPatchMap p2g(this);
2437 
2438  Array<Connection> bel_dof_list;
2439  bel_to_patch.SetSize(NumOfActiveBdrElems);
2440  bel_to_IJK.SetSize(NumOfActiveBdrElems, 2);
2441 
2442  for (int b = 0; b < GetNBP(); b++)
2443  {
2444  p2g.SetBdrPatchDofMap(b, kv, okv);
2445  const int nx = p2g.nx(); // NCP0-1
2446  const int ny = p2g.ny(); // NCP1-1
2447 
2448  // Load dofs
2449  const int nks0 = kv[0]->GetNKS();
2450  const int ord0 = kv[0]->GetOrder();
2451  const int nks1 = kv[1]->GetNKS();
2452  const int ord1 = kv[1]->GetOrder();
2453  for (int j = 0; j < nks1; j++)
2454  {
2455  if (kv[1]->isElement(j))
2456  {
2457  for (int i = 0; i < nks0; i++)
2458  {
2459  if (kv[0]->isElement(i))
2460  {
2461  if (activeBdrElem[gbe])
2462  {
2463  Connection conn(lbe,0);
2464  for (int jj = 0; jj <= ord1; jj++)
2465  {
2466  const int _jj = (okv[1] >= 0) ? (j+jj) : (ny-j-jj);
2467  for (int ii = 0; ii <= ord0; ii++)
2468  {
2469  const int _ii = (okv[0] >= 0) ? (i+ii) : (nx-i-ii);
2470  conn.to = p2g(_ii, _jj);
2471  bel_dof_list.Append(conn);
2472  }
2473  }
2474  bel_to_patch[lbe] = b;
2475  bel_to_IJK(lbe,0) = (okv[0] >= 0) ? i : (-1-i);
2476  bel_to_IJK(lbe,1) = (okv[1] >= 0) ? j : (-1-j);
2477  lbe++;
2478  }
2479  gbe++;
2480  }
2481  }
2482  }
2483  }
2484  }
2485  // We must NOT sort bel_dof_list in this case.
2486  bel_dof = new Table(NumOfActiveBdrElems, bel_dof_list);
2487 }
2488 
2489 void NURBSExtension::GetVertexLocalToGlobal(Array<int> &lvert_vert)
2490 {
2491  lvert_vert.SetSize(GetNV());
2492  for (int gv = 0; gv < GetGNV(); gv++)
2493  if (activeVert[gv] >= 0)
2494  {
2495  lvert_vert[activeVert[gv]] = gv;
2496  }
2497 }
2498 
2499 void NURBSExtension::GetElementLocalToGlobal(Array<int> &lelem_elem)
2500 {
2501  lelem_elem.SetSize(GetNE());
2502  for (int le = 0, ge = 0; ge < GetGNE(); ge++)
2503  if (activeElem[ge])
2504  {
2505  lelem_elem[le++] = ge;
2506  }
2507 }
2508 
2509 void NURBSExtension::LoadFE(int i, const FiniteElement *FE) const
2510 {
2511  const NURBSFiniteElement *NURBSFE =
2512  dynamic_cast<const NURBSFiniteElement *>(FE);
2513 
2514  if (NURBSFE->GetElement() != i)
2515  {
2516  Array<int> dofs;
2517  NURBSFE->SetIJK(el_to_IJK.GetRow(i));
2518  if (el_to_patch[i] != NURBSFE->GetPatch())
2519  {
2520  GetPatchKnotVectors(el_to_patch[i], NURBSFE->KnotVectors());
2521  NURBSFE->SetPatch(el_to_patch[i]);
2522  NURBSFE->SetOrder();
2523  }
2524  el_dof->GetRow(i, dofs);
2525  weights.GetSubVector(dofs, NURBSFE->Weights());
2526  NURBSFE->SetElement(i);
2527  }
2528 }
2529 
2530 void NURBSExtension::LoadBE(int i, const FiniteElement *BE) const
2531 {
2532  const NURBSFiniteElement *NURBSFE =
2533  dynamic_cast<const NURBSFiniteElement *>(BE);
2534 
2535  if (NURBSFE->GetElement() != i)
2536  {
2537  Array<int> dofs;
2538  NURBSFE->SetIJK(bel_to_IJK.GetRow(i));
2539  if (bel_to_patch[i] != NURBSFE->GetPatch())
2540  {
2541  GetBdrPatchKnotVectors(bel_to_patch[i], NURBSFE->KnotVectors());
2542  NURBSFE->SetPatch(bel_to_patch[i]);
2543  NURBSFE->SetOrder();
2544  }
2545  bel_dof->GetRow(i, dofs);
2546  weights.GetSubVector(dofs, NURBSFE->Weights());
2547  NURBSFE->SetElement(i);
2548  }
2549 }
2550 
2551 void NURBSExtension::ConvertToPatches(const Vector &Nodes)
2552 {
2553  delete el_dof;
2554  delete bel_dof;
2555 
2556  if (patches.Size() == 0)
2557  {
2558  GetPatchNets(Nodes, Dimension());
2559  }
2560 }
2561 
2562 void NURBSExtension::SetCoordsFromPatches(Vector &Nodes)
2563 {
2564  if (patches.Size() == 0) { return; }
2565 
2566  SetSolutionVector(Nodes, Dimension());
2567  patches.SetSize(0);
2568 }
2569 
2570 void NURBSExtension::SetKnotsFromPatches()
2571 {
2572  if (patches.Size() == 0)
2573  {
2574  mfem_error("NURBSExtension::SetKnotsFromPatches :"
2575  " No patches available!");
2576  }
2577 
2578  Array<KnotVector *> kv;
2579 
2580  for (int p = 0; p < patches.Size(); p++)
2581  {
2582  GetPatchKnotVectors(p, kv);
2583 
2584  for (int i = 0; i < kv.Size(); i++)
2585  {
2586  *kv[i] = *patches[p]->GetKV(i);
2587  }
2588  }
2589 
2590  SetOrdersFromKnotVectors();
2591 
2592  GenerateOffsets();
2593  CountElements();
2594  CountBdrElements();
2595 
2596  // all elements must be active
2597  NumOfActiveElems = NumOfElements;
2598  activeElem.SetSize(NumOfElements);
2599  activeElem = true;
2600 
2601  GenerateActiveVertices();
2602  GenerateElementDofTable();
2603  GenerateActiveBdrElems();
2604  GenerateBdrElementDofTable();
2605 }
2606 
2607 void NURBSExtension::LoadSolution(std::istream &input, GridFunction &sol) const
2608 {
2609  const FiniteElementSpace *fes = sol.FESpace();
2610  MFEM_VERIFY(fes->GetNURBSext() == this, "");
2611 
2612  sol.SetSize(fes->GetVSize());
2613 
2614  Array<const KnotVector *> kv(Dimension());
2615  NURBSPatchMap p2g(this);
2616  const int vdim = fes->GetVDim();
2617 
2618  for (int p = 0; p < GetNP(); p++)
2619  {
2620  skip_comment_lines(input, '#');
2621 
2622  p2g.SetPatchDofMap(p, kv);
2623  const int nx = kv[0]->GetNCP();
2624  const int ny = kv[1]->GetNCP();
2625  const int nz = (kv.Size() == 2) ? 1 : kv[2]->GetNCP();
2626  for (int k = 0; k < nz; k++)
2627  {
2628  for (int j = 0; j < ny; j++)
2629  {
2630  for (int i = 0; i < nx; i++)
2631  {
2632  const int l = (kv.Size() == 2) ? p2g(i,j) : p2g(i,j,k);
2633  for (int vd = 0; vd < vdim; vd++)
2634  {
2635  input >> sol(fes->DofToVDof(l,vd));
2636  }
2637  }
2638  }
2639  }
2640  }
2641 }
2642 
2643 void NURBSExtension::PrintSolution(const GridFunction &sol, std::ostream &out)
2644 const
2645 {
2646  const FiniteElementSpace *fes = sol.FESpace();
2647  MFEM_VERIFY(fes->GetNURBSext() == this, "");
2648 
2649  Array<const KnotVector *> kv(Dimension());
2650  NURBSPatchMap p2g(this);
2651  const int vdim = fes->GetVDim();
2652 
2653  for (int p = 0; p < GetNP(); p++)
2654  {
2655  out << "\n# patch " << p << "\n\n";
2656 
2657  p2g.SetPatchDofMap(p, kv);
2658  const int nx = kv[0]->GetNCP();
2659  const int ny = kv[1]->GetNCP();
2660  const int nz = (kv.Size() == 2) ? 1 : kv[2]->GetNCP();
2661  for (int k = 0; k < nz; k++)
2662  {
2663  for (int j = 0; j < ny; j++)
2664  {
2665  for (int i = 0; i < nx; i++)
2666  {
2667  const int l = (kv.Size() == 2) ? p2g(i,j) : p2g(i,j,k);
2668  out << sol(fes->DofToVDof(l,0));
2669  for (int vd = 1; vd < vdim; vd++)
2670  {
2671  out << ' ' << sol(fes->DofToVDof(l,vd));
2672  }
2673  out << '\n';
2674  }
2675  }
2676  }
2677  }
2678 }
2679 
2680 void NURBSExtension::DegreeElevate(int rel_degree, int degree)
2681 {
2682  for (int p = 0; p < patches.Size(); p++)
2683  {
2684  for (int dir = 0; dir < patches[p]->GetNKV(); dir++)
2685  {
2686  int oldd = patches[p]->GetKV(dir)->GetOrder();
2687  int newd = std::min(oldd + rel_degree, degree);
2688  if (newd > oldd)
2689  {
2690  patches[p]->DegreeElevate(dir, newd - oldd);
2691  }
2692  }
2693  }
2694 }
2695 
2696 void NURBSExtension::UniformRefinement()
2697 {
2698  for (int p = 0; p < patches.Size(); p++)
2699  {
2700  patches[p]->UniformRefinement();
2701  }
2702 }
2703 
2704 void NURBSExtension::KnotInsert(Array<KnotVector *> &kv)
2705 {
2706  Array<int> edges;
2707  Array<int> orient;
2708 
2709  Array<KnotVector *> pkv(Dimension());
2710 
2711  for (int p = 0; p < patches.Size(); p++)
2712  {
2713  patchTopo->GetElementEdges(p, edges, orient);
2714 
2715  if (Dimension()==2)
2716  {
2717  pkv[0] = kv[KnotInd(edges[0])];
2718  pkv[1] = kv[KnotInd(edges[1])];
2719  }
2720  else
2721  {
2722  pkv[0] = kv[KnotInd(edges[0])];
2723  pkv[1] = kv[KnotInd(edges[3])];
2724  pkv[2] = kv[KnotInd(edges[8])];
2725  }
2726 
2727  patches[p]->KnotInsert(pkv);
2728  }
2729 }
2730 
2731 void NURBSExtension::GetPatchNets(const Vector &coords, int vdim)
2732 {
2733  if (Dimension() == 2)
2734  {
2735  Get2DPatchNets(coords, vdim);
2736  }
2737  else
2738  {
2739  Get3DPatchNets(coords, vdim);
2740  }
2741 }
2742 
2743 void NURBSExtension::Get2DPatchNets(const Vector &coords, int vdim)
2744 {
2745  Array<const KnotVector *> kv(2);
2746  NURBSPatchMap p2g(this);
2747 
2748  patches.SetSize(GetNP());
2749  for (int p = 0; p < GetNP(); p++)
2750  {
2751  p2g.SetPatchDofMap(p, kv);
2752  patches[p] = new NURBSPatch(kv, vdim+1);
2753  NURBSPatch &Patch = *patches[p];
2754 
2755  for (int j = 0; j < kv[1]->GetNCP(); j++)
2756  {
2757  for (int i = 0; i < kv[0]->GetNCP(); i++)
2758  {
2759  const int l = p2g(i,j);
2760  for (int d = 0; d < vdim; d++)
2761  {
2762  Patch(i,j,d) = coords(l*vdim + d)*weights(l);
2763  }
2764  Patch(i,j,vdim) = weights(l);
2765  }
2766  }
2767  }
2768 }
2769 
2770 void NURBSExtension::Get3DPatchNets(const Vector &coords, int vdim)
2771 {
2772  Array<const KnotVector *> kv(3);
2773  NURBSPatchMap p2g(this);
2774 
2775  patches.SetSize(GetNP());
2776  for (int p = 0; p < GetNP(); p++)
2777  {
2778  p2g.SetPatchDofMap(p, kv);
2779  patches[p] = new NURBSPatch(kv, vdim+1);
2780  NURBSPatch &Patch = *patches[p];
2781 
2782  for (int k = 0; k < kv[2]->GetNCP(); k++)
2783  {
2784  for (int j = 0; j < kv[1]->GetNCP(); j++)
2785  {
2786  for (int i = 0; i < kv[0]->GetNCP(); i++)
2787  {
2788  const int l = p2g(i,j,k);
2789  for (int d = 0; d < vdim; d++)
2790  {
2791  Patch(i,j,k,d) = coords(l*vdim + d)*weights(l);
2792  }
2793  Patch(i,j,k,vdim) = weights(l);
2794  }
2795  }
2796  }
2797  }
2798 }
2799 
2800 void NURBSExtension::SetSolutionVector(Vector &coords, int vdim)
2801 {
2802  if (Dimension() == 2)
2803  {
2804  Set2DSolutionVector(coords, vdim);
2805  }
2806  else
2807  {
2808  Set3DSolutionVector(coords, vdim);
2809  }
2810 }
2811 
2812 void NURBSExtension::Set2DSolutionVector(Vector &coords, int vdim)
2813 {
2814  Array<const KnotVector *> kv(2);
2815  NURBSPatchMap p2g(this);
2816 
2817  weights.SetSize(GetNDof());
2818  for (int p = 0; p < GetNP(); p++)
2819  {
2820  p2g.SetPatchDofMap(p, kv);
2821  NURBSPatch &Patch = *patches[p];
2822  MFEM_ASSERT(vdim+1 == Patch.GetNC(), "");
2823 
2824  for (int j = 0; j < kv[1]->GetNCP(); j++)
2825  {
2826  for (int i = 0; i < kv[0]->GetNCP(); i++)
2827  {
2828  const int l = p2g(i,j);
2829  for (int d = 0; d < vdim; d++)
2830  {
2831  coords(l*vdim + d) = Patch(i,j,d)/Patch(i,j,vdim);
2832  }
2833  weights(l) = Patch(i,j,vdim);
2834  }
2835  }
2836  delete patches[p];
2837  }
2838 }
2839 
2840 void NURBSExtension::Set3DSolutionVector(Vector &coords, int vdim)
2841 {
2842  Array<const KnotVector *> kv(3);
2843  NURBSPatchMap p2g(this);
2844 
2845  weights.SetSize(GetNDof());
2846  for (int p = 0; p < GetNP(); p++)
2847  {
2848  p2g.SetPatchDofMap(p, kv);
2849  NURBSPatch &Patch = *patches[p];
2850  MFEM_ASSERT(vdim+1 == Patch.GetNC(), "");
2851 
2852  for (int k = 0; k < kv[2]->GetNCP(); k++)
2853  {
2854  for (int j = 0; j < kv[1]->GetNCP(); j++)
2855  {
2856  for (int i = 0; i < kv[0]->GetNCP(); i++)
2857  {
2858  const int l = p2g(i,j,k);
2859  for (int d = 0; d < vdim; d++)
2860  {
2861  coords(l*vdim + d) = Patch(i,j,k,d)/Patch(i,j,k,vdim);
2862  }
2863  weights(l) = Patch(i,j,k,vdim);
2864  }
2865  }
2866  }
2867  delete patches[p];
2868  }
2869 }
2870 
2871 
2872 #ifdef MFEM_USE_MPI
2873 ParNURBSExtension::ParNURBSExtension(const ParNURBSExtension &orig)
2874  : NURBSExtension(orig),
2875  partitioning(orig.partitioning ? new int[orig.GetGNE()] : NULL),
2876  gtopo(orig.gtopo),
2877  ldof_group(orig.ldof_group)
2878 {
2879  // Copy the partitioning, if not NULL
2880  if (partitioning)
2881  {
2882  std::memcpy(partitioning, orig.partitioning, orig.GetGNE()*sizeof(int));
2883  }
2884 }
2885 
2886 ParNURBSExtension::ParNURBSExtension(MPI_Comm comm, NURBSExtension *parent,
2887  int *part, const Array<bool> &active_bel)
2888  : gtopo(comm)
2889 {
2890  if (parent->NumOfActiveElems < parent->NumOfElements)
2891  {
2892  // SetActive (BuildGroups?) and the way the weights are copied
2893  // do not support this case
2894  mfem_error("ParNURBSExtension::ParNURBSExtension :\n"
2895  " all elements in the parent must be active!");
2896  }
2897 
2898  patchTopo = parent->patchTopo;
2899  // steal ownership of patchTopo from the 'parent' NURBS extension
2900  if (!parent->own_topo)
2901  {
2902  mfem_error("ParNURBSExtension::ParNURBSExtension :\n"
2903  " parent does not own the patch topology!");
2904  }
2905  own_topo = 1;
2906  parent->own_topo = 0;
2907 
2908  parent->edge_to_knot.Copy(edge_to_knot);
2909 
2910  parent->GetOrders().Copy(mOrders);
2911  mOrder = parent->GetOrder();
2912 
2913  NumOfKnotVectors = parent->GetNKV();
2914  knotVectors.SetSize(NumOfKnotVectors);
2915  for (int i = 0; i < NumOfKnotVectors; i++)
2916  {
2917  knotVectors[i] = new KnotVector(*parent->GetKnotVector(i));
2918  }
2919 
2920  GenerateOffsets();
2921  CountElements();
2922  CountBdrElements();
2923 
2924  // copy 'part' to 'partitioning'
2925  partitioning = new int[GetGNE()];
2926  for (int i = 0; i < GetGNE(); i++)
2927  {
2928  partitioning[i] = part[i];
2929  }
2930  SetActive(partitioning, active_bel);
2931 
2932  GenerateActiveVertices();
2933  GenerateElementDofTable();
2934  // GenerateActiveBdrElems(); // done by SetActive for now
2935  GenerateBdrElementDofTable();
2936 
2937  Table *serial_elem_dof = parent->GetElementDofTable();
2938  BuildGroups(partitioning, *serial_elem_dof);
2939 
2940  weights.SetSize(GetNDof());
2941  // copy weights from parent
2942  for (int gel = 0, lel = 0; gel < GetGNE(); gel++)
2943  {
2944  if (activeElem[gel])
2945  {
2946  int ndofs = el_dof->RowSize(lel);
2947  int *ldofs = el_dof->GetRow(lel);
2948  int *gdofs = serial_elem_dof->GetRow(gel);
2949  for (int i = 0; i < ndofs; i++)
2950  {
2951  weights(ldofs[i]) = parent->weights(gdofs[i]);
2952  }
2953  lel++;
2954  }
2955  }
2956 }
2957 
2958 ParNURBSExtension::ParNURBSExtension(NURBSExtension *parent,
2959  const ParNURBSExtension *par_parent)
2960  : gtopo(par_parent->gtopo.GetComm())
2961 {
2962  // steal all data from parent
2963  mOrder = parent->mOrder;
2964  Swap(mOrders, parent->mOrders);
2965 
2966  patchTopo = parent->patchTopo;
2967  own_topo = parent->own_topo;
2968  parent->own_topo = 0;
2969 
2970  Swap(edge_to_knot, parent->edge_to_knot);
2971 
2972  NumOfKnotVectors = parent->NumOfKnotVectors;
2973  Swap(knotVectors, parent->knotVectors);
2974 
2975  NumOfVertices = parent->NumOfVertices;
2976  NumOfElements = parent->NumOfElements;
2977  NumOfBdrElements = parent->NumOfBdrElements;
2978  NumOfDofs = parent->NumOfDofs;
2979 
2980  Swap(v_meshOffsets, parent->v_meshOffsets);
2981  Swap(e_meshOffsets, parent->e_meshOffsets);
2982  Swap(f_meshOffsets, parent->f_meshOffsets);
2983  Swap(p_meshOffsets, parent->p_meshOffsets);
2984 
2985  Swap(v_spaceOffsets, parent->v_spaceOffsets);
2986  Swap(e_spaceOffsets, parent->e_spaceOffsets);
2987  Swap(f_spaceOffsets, parent->f_spaceOffsets);
2988  Swap(p_spaceOffsets, parent->p_spaceOffsets);
2989 
2990  NumOfActiveVertices = parent->NumOfActiveVertices;
2991  NumOfActiveElems = parent->NumOfActiveElems;
2992  NumOfActiveBdrElems = parent->NumOfActiveBdrElems;
2993  NumOfActiveDofs = parent->NumOfActiveDofs;
2994 
2995  Swap(activeVert, parent->activeVert);
2996  Swap(activeElem, parent->activeElem);
2997  Swap(activeBdrElem, parent->activeBdrElem);
2998  Swap(activeDof, parent->activeDof);
2999 
3000  el_dof = parent->el_dof;
3001  bel_dof = parent->bel_dof;
3002  parent->el_dof = parent->bel_dof = NULL;
3003 
3004  Swap(el_to_patch, parent->el_to_patch);
3005  Swap(bel_to_patch, parent->bel_to_patch);
3006  Swap(el_to_IJK, parent->el_to_IJK);
3007  Swap(bel_to_IJK, parent->bel_to_IJK);
3008 
3009  Swap(weights, parent->weights);
3010  MFEM_VERIFY(!parent->HavePatches(), "");
3011 
3012  delete parent;
3013 
3014  partitioning = NULL;
3015 
3016  MFEM_VERIFY(par_parent->partitioning,
3017  "parent ParNURBSExtension has no partitioning!");
3018 
3019  // Support for the case when 'parent' is not a local NURBSExtension, i.e.
3020  // NumOfActiveElems is not the same as in 'par_parent'. In that case, we
3021  // assume 'parent' is a global NURBSExtension, i.e. all elements are active.
3022  bool extract_weights = false;
3023  if (NumOfActiveElems != par_parent->NumOfActiveElems)
3024  {
3025  MFEM_ASSERT(NumOfActiveElems == NumOfElements, "internal error");
3026 
3027  SetActive(par_parent->partitioning, par_parent->activeBdrElem);
3028  GenerateActiveVertices();
3029  delete el_dof;
3030  el_to_patch.DeleteAll();
3031  el_to_IJK.DeleteAll();
3032  GenerateElementDofTable();
3033  // GenerateActiveBdrElems(); // done by SetActive for now
3034  delete bel_dof;
3035  bel_to_patch.DeleteAll();
3036  bel_to_IJK.DeleteAll();
3037  GenerateBdrElementDofTable();
3038  extract_weights = true;
3039  }
3040 
3041  Table *glob_elem_dof = GetGlobalElementDofTable();
3042  BuildGroups(par_parent->partitioning, *glob_elem_dof);
3043  if (extract_weights)
3044  {
3045  Vector glob_weights;
3046  Swap(weights, glob_weights);
3047  weights.SetSize(GetNDof());
3048  // Copy the local 'weights' from the 'glob_weights'.
3049  // Assumption: the local element ids follow the global ordering.
3050  for (int gel = 0, lel = 0; gel < GetGNE(); gel++)
3051  {
3052  if (activeElem[gel])
3053  {
3054  int ndofs = el_dof->RowSize(lel);
3055  int *ldofs = el_dof->GetRow(lel);
3056  int *gdofs = glob_elem_dof->GetRow(gel);
3057  for (int i = 0; i < ndofs; i++)
3058  {
3059  weights(ldofs[i]) = glob_weights(gdofs[i]);
3060  }
3061  lel++;
3062  }
3063  }
3064  }
3065  delete glob_elem_dof;
3066 }
3067 
3068 Table *ParNURBSExtension::GetGlobalElementDofTable()
3069 {
3070  if (Dimension() == 2)
3071  {
3072  return Get2DGlobalElementDofTable();
3073  }
3074  else
3075  {
3076  return Get3DGlobalElementDofTable();
3077  }
3078 }
3079 
3080 Table *ParNURBSExtension::Get2DGlobalElementDofTable()
3081 {
3082  int el = 0;
3083  const KnotVector *kv[2];
3084  NURBSPatchMap p2g(this);
3085  Array<Connection> gel_dof_list;
3086 
3087  for (int p = 0; p < GetNP(); p++)
3088  {
3089  p2g.SetPatchDofMap(p, kv);
3090 
3091  // Load dofs
3092  const int ord0 = kv[0]->GetOrder();
3093  const int ord1 = kv[1]->GetOrder();
3094  for (int j = 0; j < kv[1]->GetNKS(); j++)
3095  {
3096  if (kv[1]->isElement(j))
3097  {
3098  for (int i = 0; i < kv[0]->GetNKS(); i++)
3099  {
3100  if (kv[0]->isElement(i))
3101  {
3102  Connection conn(el,0);
3103  for (int jj = 0; jj <= ord1; jj++)
3104  {
3105  for (int ii = 0; ii <= ord0; ii++)
3106  {
3107  conn.to = p2g(i+ii,j+jj);
3108  gel_dof_list.Append(conn);
3109  }
3110  }
3111  el++;
3112  }
3113  }
3114  }
3115  }
3116  }
3117  // We must NOT sort gel_dof_list in this case.
3118  return (new Table(GetGNE(), gel_dof_list));
3119 }
3120 
3121 Table *ParNURBSExtension::Get3DGlobalElementDofTable()
3122 {
3123  int el = 0;
3124  const KnotVector *kv[3];
3125  NURBSPatchMap p2g(this);
3126  Array<Connection> gel_dof_list;
3127 
3128  for (int p = 0; p < GetNP(); p++)
3129  {
3130  p2g.SetPatchDofMap(p, kv);
3131 
3132  // Load dofs
3133  const int ord0 = kv[0]->GetOrder();
3134  const int ord1 = kv[1]->GetOrder();
3135  const int ord2 = kv[2]->GetOrder();
3136  for (int k = 0; k < kv[2]->GetNKS(); k++)
3137  {
3138  if (kv[2]->isElement(k))
3139  {
3140  for (int j = 0; j < kv[1]->GetNKS(); j++)
3141  {
3142  if (kv[1]->isElement(j))
3143  {
3144  for (int i = 0; i < kv[0]->GetNKS(); i++)
3145  {
3146  if (kv[0]->isElement(i))
3147  {
3148  Connection conn(el,0);
3149  for (int kk = 0; kk <= ord2; kk++)
3150  {
3151  for (int jj = 0; jj <= ord1; jj++)
3152  {
3153  for (int ii = 0; ii <= ord0; ii++)
3154  {
3155  conn.to = p2g(i+ii,j+jj,k+kk);
3156  gel_dof_list.Append(conn);
3157  }
3158  }
3159  }
3160  el++;
3161  }
3162  }
3163  }
3164  }
3165  }
3166  }
3167  }
3168  // We must NOT sort gel_dof_list in this case.
3169  return (new Table(GetGNE(), gel_dof_list));
3170 }
3171 
3172 void ParNURBSExtension::SetActive(const int *_partitioning,
3173  const Array<bool> &active_bel)
3174 {
3175  activeElem.SetSize(GetGNE());
3176  activeElem = false;
3177  NumOfActiveElems = 0;
3178  const int MyRank = gtopo.MyRank();
3179  for (int i = 0; i < GetGNE(); i++)
3180  if (_partitioning[i] == MyRank)
3181  {
3182  activeElem[i] = true;
3183  NumOfActiveElems++;
3184  }
3185 
3186  active_bel.Copy(activeBdrElem);
3187  NumOfActiveBdrElems = 0;
3188  for (int i = 0; i < GetGNBE(); i++)
3189  if (activeBdrElem[i])
3190  {
3191  NumOfActiveBdrElems++;
3192  }
3193 }
3194 
3195 void ParNURBSExtension::BuildGroups(const int *_partitioning,
3196  const Table &elem_dof)
3197 {
3198  Table dof_proc;
3199 
3200  ListOfIntegerSets groups;
3201  IntegerSet group;
3202 
3203  Transpose(elem_dof, dof_proc); // dof_proc is dof_elem
3204  // convert elements to processors
3205  for (int i = 0; i < dof_proc.Size_of_connections(); i++)
3206  {
3207  dof_proc.GetJ()[i] = _partitioning[dof_proc.GetJ()[i]];
3208  }
3209 
3210  // the first group is the local one
3211  int MyRank = gtopo.MyRank();
3212  group.Recreate(1, &MyRank);
3213  groups.Insert(group);
3214 
3215  int dof = 0;
3216  ldof_group.SetSize(GetNDof());
3217  for (int d = 0; d < GetNTotalDof(); d++)
3218  if (activeDof[d])
3219  {
3220  group.Recreate(dof_proc.RowSize(d), dof_proc.GetRow(d));
3221  ldof_group[dof] = groups.Insert(group);
3222 
3223  dof++;
3224  }
3225 
3226  gtopo.Create(groups, 1822);
3227 }
3228 #endif // MFEM_USE_MPI
3229 
3230 
3231 void NURBSPatchMap::GetPatchKnotVectors(int p, const KnotVector *kv[])
3232 {
3233  Ext->patchTopo->GetElementVertices(p, verts);
3234  Ext->patchTopo->GetElementEdges(p, edges, oedge);
3235  if (Ext->Dimension() == 2)
3236  {
3237  kv[0] = Ext->KnotVec(edges[0]);
3238  kv[1] = Ext->KnotVec(edges[1]);
3239  }
3240  else
3241  {
3242  Ext->patchTopo->GetElementFaces(p, faces, oface);
3243 
3244  kv[0] = Ext->KnotVec(edges[0]);
3245  kv[1] = Ext->KnotVec(edges[3]);
3246  kv[2] = Ext->KnotVec(edges[8]);
3247  }
3248  opatch = 0;
3249 }
3250 
3251 void NURBSPatchMap::GetBdrPatchKnotVectors(int p, const KnotVector *kv[],
3252  int *okv)
3253 {
3254  Ext->patchTopo->GetBdrElementVertices(p, verts);
3255  Ext->patchTopo->GetBdrElementEdges(p, edges, oedge);
3256  kv[0] = Ext->KnotVec(edges[0], oedge[0], &okv[0]);
3257  if (Ext->Dimension() == 3)
3258  {
3259  faces.SetSize(1);
3260  Ext->patchTopo->GetBdrElementFace(p, &faces[0], &opatch);
3261  kv[1] = Ext->KnotVec(edges[1], oedge[1], &okv[1]);
3262  }
3263  else
3264  {
3265  opatch = oedge[0];
3266  }
3267 }
3268 
3269 void NURBSPatchMap::SetPatchVertexMap(int p, const KnotVector *kv[])
3270 {
3271  GetPatchKnotVectors(p, kv);
3272 
3273  I = kv[0]->GetNE() - 1;
3274  J = kv[1]->GetNE() - 1;
3275 
3276  for (int i = 0; i < verts.Size(); i++)
3277  {
3278  verts[i] = Ext->v_meshOffsets[verts[i]];
3279  }
3280 
3281  for (int i = 0; i < edges.Size(); i++)
3282  {
3283  edges[i] = Ext->e_meshOffsets[edges[i]];
3284  }
3285 
3286  if (Ext->Dimension() == 3)
3287  {
3288  K = kv[2]->GetNE() - 1;
3289 
3290  for (int i = 0; i < faces.Size(); i++)
3291  {
3292  faces[i] = Ext->f_meshOffsets[faces[i]];
3293  }
3294  }
3295 
3296  pOffset = Ext->p_meshOffsets[p];
3297 }
3298 
3299 void NURBSPatchMap::SetPatchDofMap(int p, const KnotVector *kv[])
3300 {
3301  GetPatchKnotVectors(p, kv);
3302 
3303  I = kv[0]->GetNCP() - 2;
3304  J = kv[1]->GetNCP() - 2;
3305 
3306  for (int i = 0; i < verts.Size(); i++)
3307  {
3308  verts[i] = Ext->v_spaceOffsets[verts[i]];
3309  }
3310 
3311  for (int i = 0; i < edges.Size(); i++)
3312  {
3313  edges[i] = Ext->e_spaceOffsets[edges[i]];
3314  }
3315 
3316  if (Ext->Dimension() == 3)
3317  {
3318  K = kv[2]->GetNCP() - 2;
3319 
3320  for (int i = 0; i < faces.Size(); i++)
3321  {
3322  faces[i] = Ext->f_spaceOffsets[faces[i]];
3323  }
3324  }
3325 
3326  pOffset = Ext->p_spaceOffsets[p];
3327 }
3328 
3329 void NURBSPatchMap::SetBdrPatchVertexMap(int p, const KnotVector *kv[],
3330  int *okv)
3331 {
3332  GetBdrPatchKnotVectors(p, kv, okv);
3333 
3334  I = kv[0]->GetNE() - 1;
3335 
3336  for (int i = 0; i < verts.Size(); i++)
3337  {
3338  verts[i] = Ext->v_meshOffsets[verts[i]];
3339  }
3340 
3341  if (Ext->Dimension() == 2)
3342  {
3343  pOffset = Ext->e_meshOffsets[edges[0]];
3344  }
3345  else
3346  {
3347  J = kv[1]->GetNE() - 1;
3348 
3349  for (int i = 0; i < edges.Size(); i++)
3350  {
3351  edges[i] = Ext->e_meshOffsets[edges[i]];
3352  }
3353 
3354  pOffset = Ext->f_meshOffsets[faces[0]];
3355  }
3356 }
3357 
3358 void NURBSPatchMap::SetBdrPatchDofMap(int p, const KnotVector *kv[], int *okv)
3359 {
3360  GetBdrPatchKnotVectors(p, kv, okv);
3361 
3362  I = kv[0]->GetNCP() - 2;
3363 
3364  for (int i = 0; i < verts.Size(); i++)
3365  {
3366  verts[i] = Ext->v_spaceOffsets[verts[i]];
3367  }
3368 
3369  if (Ext->Dimension() == 2)
3370  {
3371  pOffset = Ext->e_spaceOffsets[edges[0]];
3372  }
3373  else
3374  {
3375  J = kv[1]->GetNCP() - 2;
3376 
3377  for (int i = 0; i < edges.Size(); i++)
3378  {
3379  edges[i] = Ext->e_spaceOffsets[edges[i]];
3380  }
3381 
3382  pOffset = Ext->f_spaceOffsets[faces[0]];
3383  }
3384 }
3385 
3386 }
KnotVector * DegreeElevate(int t) const
Definition: nurbs.cpp:57
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 KnotInsert(int dir, const KnotVector &knot)
Definition: nurbs.cpp:559
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:328
void Rotate3D(double normal[], double angle)
Definition: nurbs.cpp:1024
void GetElements()
Count the number of elements.
Definition: nurbs.cpp:100
int NumOfControlPoints
Definition: nurbs.hpp:37
void init(int dim_)
Definition: nurbs.cpp:277
Vector knot
Definition: nurbs.hpp:36
void SwapDirections(int dir1, int dir2)
Definition: nurbs.cpp:960
KnotVector()
Create KnotVector.
Definition: nurbs.hpp:41
static const int MaxOrder
Definition: nurbs.hpp:34
void UniformRefinement(Vector &newknots) const
Definition: nurbs.cpp:86
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
KnotVector & operator=(const KnotVector &kv)
Definition: nurbs.cpp:45
void CalcShape(Vector &shape, int i, double xi) const
Definition: nurbs.cpp:132
NURBSPatch * Revolve3D(NURBSPatch &patch, double n[], double ang, int times)
Definition: nurbs.cpp:1120
int GetNCP() const
Definition: nurbs.hpp:50
void Difference(const KnotVector &kv, Vector &diff) const
Definition: nurbs.cpp:243
Vector data type.
Definition: vector.hpp:48
void Print(std::ostream &out) const
Definition: nurbs.cpp:125
void CalcDShape(Vector &grad, int i, double xi) const
Definition: nurbs.cpp:158
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
int SetLoopDirection(int dir)
Definition: nurbs.cpp:482