MFEM  v4.6.0 Finite element discretization library
nurbs.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11
13 #include "../fem/fem.hpp"
14 #include "../general/text.hpp"
15
16 #include <fstream>
17 #include <algorithm>
18 #if defined(_MSC_VER) && (_MSC_VER < 1800)
19 #include <float.h>
20 #define copysign _copysign
21 #endif
22
23 namespace mfem
24 {
25
26 using namespace std;
27
28 const int KnotVector::MaxOrder = 10;
29
30 KnotVector::KnotVector(std::istream &input)
31 {
32  input >> Order >> NumOfControlPoints;
33  knot.Load(input, NumOfControlPoints + Order + 1);
34  GetElements();
35 }
36
37 KnotVector::KnotVector(int Order_, int NCP)
38 {
39  Order = Order_;
40  NumOfControlPoints = NCP;
41  knot.SetSize(NumOfControlPoints + Order + 1);
42
43  knot = -1.;
44 }
45
47 {
48  Order = kv.Order;
49  NumOfControlPoints = kv.NumOfControlPoints;
50  NumOfElements = kv.NumOfElements;
51  knot = kv.knot;
52  // alternatively, re-compute NumOfElements
53  // GetElements();
54
55  return *this;
56 }
57
59 {
60  if (t < 0)
61  {
62  mfem_error("KnotVector::DegreeElevate :\n"
63  " Parent KnotVector order higher than child");
64  }
65
66  const int nOrder = Order + t;
67  KnotVector *newkv = new KnotVector(nOrder, GetNCP() + t);
68
69  for (int i = 0; i <= nOrder; i++)
70  {
71  (*newkv)[i] = knot(0);
72  }
73  for (int i = nOrder + 1; i < newkv->GetNCP(); i++)
74  {
75  (*newkv)[i] = knot(i - t);
76  }
77  for (int i = 0; i <= nOrder; i++)
78  {
79  (*newkv)[newkv->GetNCP() + i] = knot(knot.Size()-1);
80  }
81
82  newkv->GetElements();
83
84  return newkv;
85 }
86
88 {
89  newknots.SetSize(NumOfElements);
90  int j = 0;
91  for (int i = 0; i < knot.Size()-1; i++)
92  {
93  if (knot(i) != knot(i+1))
94  {
95  newknots(j) = 0.5*(knot(i) + knot(i+1));
96  j++;
97  }
98  }
99 }
100
102 {
103  NumOfElements = 0;
104  for (int i = Order; i < NumOfControlPoints; i++)
105  {
106  if (knot(i) != knot(i+1))
107  {
108  NumOfElements++;
109  }
110  }
111 }
112
114 {
115  double apb = knot(0) + knot(knot.Size()-1);
116
117  int ns = (NumOfControlPoints - Order)/2;
118  for (int i = 1; i <= ns; i++)
119  {
120  double tmp = apb - knot(Order + i);
121  knot(Order + i) = apb - knot(NumOfControlPoints - i);
122  knot(NumOfControlPoints - i) = tmp;
123  }
124 }
125
126 void KnotVector::Print(std::ostream &os) const
127 {
128  os << Order << ' ' << NumOfControlPoints << ' ';
129  knot.Print(os, knot.Size());
130 }
131
132
133 void KnotVector::PrintFunctions(std::ostream &os, int samples) const
134 {
135  Vector shape(Order+1);
136
137  double x, dx = 1.0/double (samples - 1);
138
139  for (int i = 0; i <GetNE() ; i++)
140  {
141  for (int j = 0; j <samples; j++)
142  {
143  x =j*dx;
144  os<< x + i;
145
146  CalcShape ( shape, i, x);
147  for (int d = 0; d < Order+1; d++) { os<<"\t"<<shape[d]; }
148
149  CalcDShape ( shape, i, x);
150  for (int d = 0; d < Order+1; d++) { os<<"\t"<<shape[d]; }
151
152  CalcD2Shape ( shape, i, x);
153  for (int d = 0; d < Order+1; d++) { os<<"\t"<<shape[d]; }
154  os<<endl;
155  }
156  }
157 }
158
159 // Routine from "The NURBS book" - 2nd ed - Piegl and Tiller
160 // Algorithm A2.2 p. 70
161 void KnotVector::CalcShape(Vector &shape, int i, double xi) const
162 {
163  MFEM_ASSERT(Order <= MaxOrder, "Order > MaxOrder!");
164
165  int p = Order;
166  int ip = (i >= 0) ? (i + p) : (-1 - i + p);
167  double u = getKnotLocation((i >= 0) ? xi : 1. - xi, ip), saved, tmp;
168  double left[MaxOrder+1], right[MaxOrder+1];
169
170  shape(0) = 1.;
171  for (int j = 1; j <= p; ++j)
172  {
173  left[j] = u - knot(ip+1-j);
174  right[j] = knot(ip+j) - u;
175  saved = 0.;
176  for (int r = 0; r < j; ++r)
177  {
178  tmp = shape(r)/(right[r+1] + left[j-r]);
179  shape(r) = saved + right[r+1]*tmp;
180  saved = left[j-r]*tmp;
181  }
182  shape(j) = saved;
183  }
184 }
185
186 // Routine from "The NURBS book" - 2nd ed - Piegl and Tiller
187 // Algorithm A2.3 p. 72
188 void KnotVector::CalcDShape(Vector &grad, int i, double xi) const
189 {
190  int p = Order, rk, pk;
191  int ip = (i >= 0) ? (i + p) : (-1 - i + p);
192  double u = getKnotLocation((i >= 0) ? xi : 1. - xi, ip), temp, saved, d;
193  double ndu[MaxOrder+1][MaxOrder+1], left[MaxOrder+1], right[MaxOrder+1];
194
195 #ifdef MFEM_DEBUG
196  if (p > MaxOrder)
197  {
198  mfem_error("KnotVector::CalcDShape : Order > MaxOrder!");
199  }
200 #endif
201
202  ndu[0][0] = 1.0;
203  for (int j = 1; j <= p; j++)
204  {
205  left[j] = u - knot(ip-j+1);
206  right[j] = knot(ip+j) - u;
207  saved = 0.0;
208  for (int r = 0; r < j; r++)
209  {
210  ndu[j][r] = right[r+1] + left[j-r];
211  temp = ndu[r][j-1]/ndu[j][r];
212  ndu[r][j] = saved + right[r+1]*temp;
213  saved = left[j-r]*temp;
214  }
215  ndu[j][j] = saved;
216  }
217
218  for (int r = 0; r <= p; ++r)
219  {
220  d = 0.0;
221  rk = r-1;
222  pk = p-1;
223  if (r >= 1)
224  {
225  d = ndu[rk][pk]/ndu[p][rk];
226  }
227  if (r <= pk)
228  {
229  d -= ndu[r][pk]/ndu[p][r];
230  }
232  }
233
234  if (i >= 0)
235  {
236  grad *= p*(knot(ip+1) - knot(ip));
237  }
238  else
239  {
240  grad *= p*(knot(ip) - knot(ip+1));
241  }
242 }
243
244 // Routine from "The NURBS book" - 2nd ed - Piegl and Tiller
245 void KnotVector::CalcDnShape(Vector &gradn, int n, int i, double xi) const
246 {
247  int p = Order, rk, pk, j1, j2,r,j,k;
248  int ip = (i >= 0) ? (i + p) : (-1 - i + p);
249  double u = getKnotLocation((i >= 0) ? xi : 1. - xi, ip);
250  double temp, saved, d;
251  double a[2][MaxOrder+1],ndu[MaxOrder+1][MaxOrder+1], left[MaxOrder+1],
252  right[MaxOrder+1];
253
254 #ifdef MFEM_DEBUG
255  if (p > MaxOrder)
256  {
257  mfem_error("KnotVector::CalcDnShape : Order > MaxOrder!");
258  }
259 #endif
260
261  ndu[0][0] = 1.0;
262  for (j = 1; j <= p; j++)
263  {
264  left[j] = u - knot(ip-j+1);
265  right[j] = knot(ip+j)- u;
266
267  saved = 0.0;
268  for (r = 0; r < j; r++)
269  {
270  ndu[j][r] = right[r+1] + left[j-r];
271  temp = ndu[r][j-1]/ndu[j][r];
272  ndu[r][j] = saved + right[r+1]*temp;
273  saved = left[j-r]*temp;
274  }
275  ndu[j][j] = saved;
276  }
277
278  for (r = 0; r <= p; r++)
279  {
280  int s1 = 0;
281  int s2 = 1;
282  a[0][0] = 1.0;
283  for (k = 1; k <= n; k++)
284  {
285  d = 0.0;
286  rk = r-k;
287  pk = p-k;
288  if (r >= k)
289  {
290  a[s2][0] = a[s1][0]/ndu[pk+1][rk];
291  d = a[s2][0]*ndu[rk][pk];
292  }
293
294  if (rk >= -1)
295  {
296  j1 = 1;
297  }
298  else
299  {
300  j1 = -rk;
301  }
302
303  if (r-1<= pk)
304  {
305  j2 = k-1;
306  }
307  else
308  {
309  j2 = p-r;
310  }
311
312  for (j = j1; j <= j2; j++)
313  {
314  a[s2][j] = (a[s1][j] - a[s1][j-1])/ndu[pk+1][rk+j];
315  d += a[s2][j]*ndu[rk+j][pk];
316  }
317
318  if (r <= pk)
319  {
320  a[s2][k] = - a[s1][k-1]/ndu[pk+1][r];
321  d += a[s2][j]*ndu[rk+j][pk];
322  }
324  j = s1;
325  s1 = s2;
326  s2 = j;
327  }
328  }
329
330  if (i >= 0)
331  {
332  u = (knot(ip+1) - knot(ip));
333  }
334  else
335  {
336  u = (knot(ip) - knot(ip+1));
337  }
338
339  temp = p*u;
340  for (k = 1; k <= n-1; k++) { temp *= (p-k)*u; }
341
342  for (j = 0; j <= p; j++) { gradn[j] *= temp; }
343
344 }
345
346 void KnotVector::FindMaxima(Array<int> &ks,
347  Vector &xi,
348  Vector &u)
349 {
350  Vector shape(Order+1);
351  Vector maxima(GetNCP());
352  double arg1, arg2, arg, max1, max2, max;
353
354  xi.SetSize(GetNCP());
355  u.SetSize(GetNCP());
356  ks.SetSize(GetNCP());
357  for (int j = 0; j <GetNCP(); j++)
358  {
359  maxima[j] = 0;
360  for (int d = 0; d < Order+1; d++)
361  {
362  int i = j - d;
363  if (isElement(i))
364  {
365  arg1 = 1e-16;
366  CalcShape(shape, i, arg1);
367  max1 = shape[d];
368
369  arg2 = 1-(1e-16);
370  CalcShape(shape, i, arg2);
371  max2 = shape[d];
372
373  arg = (arg1 + arg2)/2;
374  CalcShape(shape, i, arg);
375  max = shape[d];
376
377  while ( ( max > max1 ) || (max > max2) )
378  {
379  if (max1 < max2)
380  {
381  max1 = max;
382  arg1 = arg;
383  }
384  else
385  {
386  max2 = max;
387  arg2 = arg;
388  }
389
390  arg = (arg1 + arg2)/2;
391  CalcShape ( shape, i, arg);
392  max = shape[d];
393  }
394
395  if (max > maxima[j])
396  {
397  maxima[j] = max;
398  ks[j] = i;
399  xi[j] = arg;
400  u[j] = getKnotLocation(arg, i+Order);
401  }
402  }
403  }
404  }
405 }
406
407 // Routine from "The NURBS book" - 2nd ed - Piegl and Tiller
408 // Algorithm A9.1 p. 369
409 void KnotVector::FindInterpolant(Array<Vector*> &x)
410 {
411  int order = GetOrder();
412  int ncp = GetNCP();
413
414  // Find interpolation points
415  Vector xi_args, u_args;
416  Array<int> i_args;
417  FindMaxima(i_args,xi_args, u_args);
418
419  // Assemble collocation matrix
420  Vector shape(order+1);
421  DenseMatrix A(ncp,ncp);
422  A = 0.0;
423  for (int i = 0; i < ncp; i++)
424  {
425  CalcShape ( shape, i_args[i], xi_args[i]);
426  for (int p = 0; p < order+1; p++)
427  {
428  A(i,i_args[i] + p) = shape[p];
429  }
430  }
431
432  // Solve problems
433  A.Invert();
434  Vector tmp;
435  for (int i= 0; i < x.Size(); i++)
436  {
437  tmp = *x[i];
438  A.Mult(tmp,*x[i]);
439  }
440 }
441
442 int KnotVector::findKnotSpan(double u) const
443 {
444  int low, mid, high;
445
446  if (u == knot(NumOfControlPoints+Order))
447  {
448  mid = NumOfControlPoints;
449  }
450  else
451  {
452  low = Order;
453  high = NumOfControlPoints + 1;
454  mid = (low + high)/2;
455  while ( (u < knot(mid-1)) || (u > knot(mid)) )
456  {
457  if (u < knot(mid-1))
458  {
459  high = mid;
460  }
461  else
462  {
463  low = mid;
464  }
465  mid = (low + high)/2;
466  }
467  }
468  return mid;
469 }
470
471 void KnotVector::Difference(const KnotVector &kv, Vector &diff) const
472 {
473  if (Order != kv.GetOrder())
474  {
475  mfem_error("KnotVector::Difference :\n"
476  " Can not compare knot vectors with different orders!");
477  }
478
479  int s = kv.Size() - Size();
480  if (s < 0)
481  {
482  kv.Difference(*this, diff);
483  return;
484  }
485
486  diff.SetSize(s);
487
488  s = 0;
489  int i = 0;
490  for (int j = 0; j < kv.Size(); j++)
491  {
492  if (abs(knot(i) - kv[j]) < 2 * std::numeric_limits<double>::epsilon())
493  {
494  i++;
495  }
496  else
497  {
498  diff(s) = kv[j];
499  s++;
500  }
501  }
502 }
503
504 void NURBSPatch::init(int dim_)
505 {
506  Dim = dim_;
507  sd = nd = -1;
508
509  if (kv.Size() == 1)
510  {
511  ni = kv[0]->GetNCP();
512  nj = -1;
513  nk = -1;
514
515  data = new double[ni*Dim];
516
517 #ifdef MFEM_DEBUG
518  for (int i = 0; i < ni*Dim; i++)
519  {
520  data[i] = -999.99;
521  }
522 #endif
523  }
524  else if (kv.Size() == 2)
525  {
526  ni = kv[0]->GetNCP();
527  nj = kv[1]->GetNCP();
528  nk = -1;
529
530  data = new double[ni*nj*Dim];
531
532 #ifdef MFEM_DEBUG
533  for (int i = 0; i < ni*nj*Dim; i++)
534  {
535  data[i] = -999.99;
536  }
537 #endif
538  }
539  else if (kv.Size() == 3)
540  {
541  ni = kv[0]->GetNCP();
542  nj = kv[1]->GetNCP();
543  nk = kv[2]->GetNCP();
544
545  data = new double[ni*nj*nk*Dim];
546
547 #ifdef MFEM_DEBUG
548  for (int i = 0; i < ni*nj*nk*Dim; i++)
549  {
550  data[i] = -999.99;
551  }
552 #endif
553  }
554  else
555  {
556  mfem_error("NURBSPatch::init : Wrong dimension of knotvectors!");
557  }
558 }
559
560 NURBSPatch::NURBSPatch(const NURBSPatch &orig)
561  : ni(orig.ni), nj(orig.nj), nk(orig.nk), Dim(orig.Dim),
562  data(NULL), kv(orig.kv.Size()), nd(orig.nd), ls(orig.ls), sd(orig.sd)
563 {
564  // Allocate and copy data:
565  const int data_size = Dim*ni*nj*((kv.Size() == 2) ? 1 : nk);
566  data = new double[data_size];
567  std::memcpy(data, orig.data, data_size*sizeof(double));
568
569  // Copy the knot vectors:
570  for (int i = 0; i < kv.Size(); i++)
571  {
572  kv[i] = new KnotVector(*orig.kv[i]);
573  }
574 }
575
576 NURBSPatch::NURBSPatch(std::istream &input)
577 {
578  int pdim, dim, size = 1;
579  string ident;
580
581  input >> ws >> ident >> pdim; // knotvectors
582  kv.SetSize(pdim);
583  for (int i = 0; i < pdim; i++)
584  {
585  kv[i] = new KnotVector(input);
586  size *= kv[i]->GetNCP();
587  }
588
589  input >> ws >> ident >> dim; // dimension
590  init(dim + 1);
591
592  input >> ws >> ident; // controlpoints (homogeneous coordinates)
593  if (ident == "controlpoints" || ident == "controlpoints_homogeneous")
594  {
595  for (int j = 0, i = 0; i < size; i++)
596  for (int d = 0; d <= dim; d++, j++)
597  {
598  input >> data[j];
599  }
600  }
601  else // "controlpoints_cartesian" (Cartesian coordinates with weight)
602  {
603  for (int j = 0, i = 0; i < size; i++)
604  {
605  for (int d = 0; d <= dim; d++)
606  {
607  input >> data[j+d];
608  }
609  for (int d = 0; d < dim; d++)
610  {
611  data[j+d] *= data[j+dim];
612  }
613  j += (dim+1);
614  }
615  }
616 }
617
618 NURBSPatch::NURBSPatch(const KnotVector *kv0, const KnotVector *kv1, int dim_)
619 {
620  kv.SetSize(2);
621  kv[0] = new KnotVector(*kv0);
622  kv[1] = new KnotVector(*kv1);
623  init(dim_);
624 }
625
626 NURBSPatch::NURBSPatch(const KnotVector *kv0, const KnotVector *kv1,
627  const KnotVector *kv2, int dim_)
628 {
629  kv.SetSize(3);
630  kv[0] = new KnotVector(*kv0);
631  kv[1] = new KnotVector(*kv1);
632  kv[2] = new KnotVector(*kv2);
633  init(dim_);
634 }
635
636 NURBSPatch::NURBSPatch(Array<const KnotVector *> &kv_, int dim_)
637 {
638  kv.SetSize(kv_.Size());
639  for (int i = 0; i < kv.Size(); i++)
640  {
641  kv[i] = new KnotVector(*kv_[i]);
642  }
643  init(dim_);
644 }
645
646 NURBSPatch::NURBSPatch(NURBSPatch *parent, int dir, int Order, int NCP)
647 {
648  kv.SetSize(parent->kv.Size());
649  for (int i = 0; i < kv.Size(); i++)
650  if (i != dir)
651  {
652  kv[i] = new KnotVector(*parent->kv[i]);
653  }
654  else
655  {
656  kv[i] = new KnotVector(Order, NCP);
657  }
658  init(parent->Dim);
659 }
660
661 void NURBSPatch::swap(NURBSPatch *np)
662 {
663  if (data != NULL)
664  {
665  delete [] data;
666  }
667
668  for (int i = 0; i < kv.Size(); i++)
669  {
670  if (kv[i]) { delete kv[i]; }
671  }
672
673  data = np->data;
674  np->kv.Copy(kv);
675
676  ni = np->ni;
677  nj = np->nj;
678  nk = np->nk;
679  Dim = np->Dim;
680
681  np->data = NULL;
682  np->kv.SetSize(0);
683
684  delete np;
685 }
686
687 NURBSPatch::~NURBSPatch()
688 {
689  if (data != NULL)
690  {
691  delete [] data;
692  }
693
694  for (int i = 0; i < kv.Size(); i++)
695  {
696  if (kv[i]) { delete kv[i]; }
697  }
698 }
699
700 void NURBSPatch::Print(std::ostream &os) const
701 {
702  int size = 1;
703
704  os << "knotvectors\n" << kv.Size() << '\n';
705  for (int i = 0; i < kv.Size(); i++)
706  {
707  kv[i]->Print(os);
708  size *= kv[i]->GetNCP();
709  }
710
711  os << "\ndimension\n" << Dim - 1
712  << "\n\ncontrolpoints\n";
713  for (int j = 0, i = 0; i < size; i++)
714  {
715  os << data[j++];
716  for (int d = 1; d < Dim; d++)
717  {
718  os << ' ' << data[j++];
719  }
720  os << '\n';
721  }
722 }
723
724 int NURBSPatch::SetLoopDirection(int dir)
725 {
726  if (nj == -1)
727  {
728  if (dir == 0)
729  {
730  sd = Dim;
731  nd = ni;
732  ls = Dim;
733  return ls;
734  }
735  else
736  {
737  mfem::err << "NURBSPatch::SetLoopDirection :\n"
738  " Direction error in 1D patch, dir = " << dir << '\n';
739  mfem_error();
740  }
741  }
742  else if (nk == -1)
743  {
744  if (dir == 0)
745  {
746  sd = Dim;
747  nd = ni;
748  ls = nj*Dim;
749  return ls;
750  }
751  else if (dir == 1)
752  {
753  sd = ni*Dim;
754  nd = nj;
755  ls = ni*Dim;
756  return ls;
757  }
758  else
759  {
760  mfem::err << "NURBSPatch::SetLoopDirection :\n"
761  " Direction error in 2D patch, dir = " << dir << '\n';
762  mfem_error();
763  }
764  }
765  else
766  {
767  if (dir == 0)
768  {
769  sd = Dim;
770  nd = ni;
771  ls = nj*nk*Dim;
772  return ls;
773  }
774  else if (dir == 1)
775  {
776  sd = ni*Dim;
777  nd = nj;
778  ls = ni*nk*Dim;
779  return ls;
780  }
781  else if (dir == 2)
782  {
783  sd = ni*nj*Dim;
784  nd = nk;
785  ls = ni*nj*Dim;
786  return ls;
787  }
788  else
789  {
790  mfem::err << "NURBSPatch::SetLoopDirection :\n"
791  " Direction error in 3D patch, dir = " << dir << '\n';
792  mfem_error();
793  }
794  }
795
796  return -1;
797 }
798
799 void NURBSPatch::UniformRefinement()
800 {
801  Vector newknots;
802  for (int dir = 0; dir < kv.Size(); dir++)
803  {
804  kv[dir]->UniformRefinement(newknots);
805  KnotInsert(dir, newknots);
806  }
807 }
808
809 void NURBSPatch::KnotInsert(Array<KnotVector *> &newkv)
810 {
811  for (int dir = 0; dir < kv.Size(); dir++)
812  {
813  KnotInsert(dir, *newkv[dir]);
814  }
815 }
816
817 void NURBSPatch::KnotInsert(int dir, const KnotVector &newkv)
818 {
819  if (dir >= kv.Size() || dir < 0)
820  {
821  mfem_error("NURBSPatch::KnotInsert : Incorrect direction!");
822  }
823
824  int t = newkv.GetOrder() - kv[dir]->GetOrder();
825
826  if (t > 0)
827  {
828  DegreeElevate(dir, t);
829  }
830  else if (t < 0)
831  {
832  mfem_error("NURBSPatch::KnotInsert : Incorrect order!");
833  }
834
835  Vector diff;
836  GetKV(dir)->Difference(newkv, diff);
837  if (diff.Size() > 0)
838  {
839  KnotInsert(dir, diff);
840  }
841 }
842
843 void NURBSPatch::KnotInsert(Array<Vector *> &newkv)
844 {
845  for (int dir = 0; dir < kv.Size(); dir++)
846  {
847  KnotInsert(dir, *newkv[dir]);
848  }
849 }
850
851 // Routine from "The NURBS book" - 2nd ed - Piegl and Tiller
852 void NURBSPatch::KnotInsert(int dir, const Vector &knot)
853 {
854  if (knot.Size() == 0 ) { return; }
855
856  if (dir >= kv.Size() || dir < 0)
857  {
858  mfem_error("NURBSPatch::KnotInsert : Incorrect direction!");
859  }
860
861  NURBSPatch &oldp = *this;
862  KnotVector &oldkv = *kv[dir];
863
864  NURBSPatch *newpatch = new NURBSPatch(this, dir, oldkv.GetOrder(),
865  oldkv.GetNCP() + knot.Size());
866  NURBSPatch &newp = *newpatch;
867  KnotVector &newkv = *newp.GetKV(dir);
868
869  int size = oldp.SetLoopDirection(dir);
870  if (size != newp.SetLoopDirection(dir))
871  {
872  mfem_error("NURBSPatch::KnotInsert : Size mismatch!");
873  }
874
875  int rr = knot.Size() - 1;
876  int a = oldkv.findKnotSpan(knot(0)) - 1;
877  int b = oldkv.findKnotSpan(knot(rr)) - 1;
878  int pl = oldkv.GetOrder();
879  int ml = oldkv.GetNCP();
880
881  for (int j = 0; j <= a; j++)
882  {
883  newkv[j] = oldkv[j];
884  }
885  for (int j = b+pl; j <= ml+pl; j++)
886  {
887  newkv[j+rr+1] = oldkv[j];
888  }
889  for (int k = 0; k <= (a-pl); k++)
890  {
891  for (int ll = 0; ll < size; ll++)
892  {
893  newp.slice(k,ll) = oldp.slice(k,ll);
894  }
895  }
896  for (int k = (b-1); k < ml; k++)
897  {
898  for (int ll = 0; ll < size; ll++)
899  {
900  newp.slice(k+rr+1,ll) = oldp.slice(k,ll);
901  }
902  }
903
904  int i = b+pl-1;
905  int k = b+pl+rr;
906
907  for (int j = rr; j >= 0; j--)
908  {
909  while ( (knot(j) <= oldkv[i]) && (i > a) )
910  {
911  newkv[k] = oldkv[i];
912  for (int ll = 0; ll < size; ll++)
913  {
914  newp.slice(k-pl-1,ll) = oldp.slice(i-pl-1,ll);
915  }
916
917  k--;
918  i--;
919  }
920
921  for (int ll = 0; ll < size; ll++)
922  {
923  newp.slice(k-pl-1,ll) = newp.slice(k-pl,ll);
924  }
925
926  for (int l = 1; l <= pl; l++)
927  {
928  int ind = k-pl+l;
929  double alfa = newkv[k+l] - knot(j);
930  if (fabs(alfa) == 0.0)
931  {
932  for (int ll = 0; ll < size; ll++)
933  {
934  newp.slice(ind-1,ll) = newp.slice(ind,ll);
935  }
936  }
937  else
938  {
939  alfa = alfa/(newkv[k+l] - oldkv[i-pl+l]);
940  for (int ll = 0; ll < size; ll++)
941  {
942  newp.slice(ind-1,ll) = alfa*newp.slice(ind-1,ll) + (1.0-alfa)*newp.slice(ind,
943  ll);
944  }
945  }
946  }
947
948  newkv[k] = knot(j);
949  k--;
950  }
951
952  newkv.GetElements();
953
954  swap(newpatch);
955 }
956
957 void NURBSPatch::DegreeElevate(int t)
958 {
959  for (int dir = 0; dir < kv.Size(); dir++)
960  {
961  DegreeElevate(dir, t);
962  }
963 }
964
965 // Routine from "The NURBS book" - 2nd ed - Piegl and Tiller
966 void NURBSPatch::DegreeElevate(int dir, int t)
967 {
968  if (dir >= kv.Size() || dir < 0)
969  {
970  mfem_error("NURBSPatch::DegreeElevate : Incorrect direction!");
971  }
972
973  int i, j, k, kj, mpi, mul, mh, kind, cind, first, last;
974  int r, a, b, oldr, save, s, tr, lbz, rbz, l;
975  double inv, ua, ub, numer, alf, den, bet, gam;
976
977  NURBSPatch &oldp = *this;
978  KnotVector &oldkv = *kv[dir];
979  oldkv.GetElements();
980
981  NURBSPatch *newpatch = new NURBSPatch(this, dir, oldkv.GetOrder() + t,
982  oldkv.GetNCP() + oldkv.GetNE()*t);
983  NURBSPatch &newp = *newpatch;
984  KnotVector &newkv = *newp.GetKV(dir);
985
986  int size = oldp.SetLoopDirection(dir);
987  if (size != newp.SetLoopDirection(dir))
988  {
989  mfem_error("NURBSPatch::DegreeElevate : Size mismatch!");
990  }
991
992  int p = oldkv.GetOrder();
993  int n = oldkv.GetNCP()-1;
994
995  DenseMatrix bezalfs (p+t+1, p+1);
996  DenseMatrix bpts (p+1, size);
997  DenseMatrix ebpts (p+t+1, size);
998  DenseMatrix nextbpts(p-1, size);
999  Vector alphas (p-1);
1000
1001  int m = n + p + 1;
1002  int ph = p + t;
1003  int ph2 = ph/2;
1004
1005  {
1006  Array2D<int> binom(ph+1, ph+1);
1007  for (i = 0; i <= ph; i++)
1008  {
1009  binom(i,0) = binom(i,i) = 1;
1010  for (j = 1; j < i; j++)
1011  {
1012  binom(i,j) = binom(i-1,j) + binom(i-1,j-1);
1013  }
1014  }
1015
1016  bezalfs(0,0) = 1.0;
1017  bezalfs(ph,p) = 1.0;
1018
1019  for (i = 1; i <= ph2; i++)
1020  {
1021  inv = 1.0/binom(ph,i);
1022  mpi = min(p,i);
1023  for (j = max(0,i-t); j <= mpi; j++)
1024  {
1025  bezalfs(i,j) = inv*binom(p,j)*binom(t,i-j);
1026  }
1027  }
1028  }
1029
1030  for (i = ph2+1; i < ph; i++)
1031  {
1032  mpi = min(p,i);
1033  for (j = max(0,i-t); j <= mpi; j++)
1034  {
1035  bezalfs(i,j) = bezalfs(ph-i,p-j);
1036  }
1037  }
1038
1039  mh = ph;
1040  kind = ph + 1;
1041  r = -1;
1042  a = p;
1043  b = p + 1;
1044  cind = 1;
1045  ua = oldkv[0];
1046  for (l = 0; l < size; l++)
1047  {
1048  newp.slice(0,l) = oldp.slice(0,l);
1049  }
1050  for (i = 0; i <= ph; i++)
1051  {
1052  newkv[i] = ua;
1053  }
1054
1055  for (i = 0; i <= p; i++)
1056  {
1057  for (l = 0; l < size; l++)
1058  {
1059  bpts(i,l) = oldp.slice(i,l);
1060  }
1061  }
1062
1063  while (b < m)
1064  {
1065  i = b;
1066  while (b < m && oldkv[b] == oldkv[b+1]) { b++; }
1067
1068  mul = b-i+1;
1069
1070  mh = mh + mul + t;
1071  ub = oldkv[b];
1072  oldr = r;
1073  r = p-mul;
1074  if (oldr > 0) { lbz = (oldr+2)/2; }
1075  else { lbz = 1; }
1076
1077  if (r > 0) { rbz = ph-(r+1)/2; }
1078  else { rbz = ph; }
1079
1080  if (r > 0)
1081  {
1082  numer = ub - ua;
1083  for (k = p ; k > mul; k--)
1084  {
1085  alphas[k-mul-1] = numer/(oldkv[a+k]-ua);
1086  }
1087
1088  for (j = 1; j <= r; j++)
1089  {
1090  save = r-j;
1091  s = mul+j;
1092  for (k = p; k >= s; k--)
1093  {
1094  for (l = 0; l < size; l++)
1095  bpts(k,l) = (alphas[k-s]*bpts(k,l) +
1096  (1.0-alphas[k-s])*bpts(k-1,l));
1097  }
1098  for (l = 0; l < size; l++)
1099  {
1100  nextbpts(save,l) = bpts(p,l);
1101  }
1102  }
1103  }
1104
1105  for (i = lbz; i <= ph; i++)
1106  {
1107  for (l = 0; l < size; l++)
1108  {
1109  ebpts(i,l) = 0.0;
1110  }
1111  mpi = min(p,i);
1112  for (j = max(0,i-t); j <= mpi; j++)
1113  {
1114  for (l = 0; l < size; l++)
1115  {
1116  ebpts(i,l) += bezalfs(i,j)*bpts(j,l);
1117  }
1118  }
1119  }
1120
1121  if (oldr > 1)
1122  {
1123  first = kind-2;
1124  last = kind;
1125  den = ub-ua;
1126  bet = (ub-newkv[kind-1])/den;
1127
1128  for (tr = 1; tr < oldr; tr++)
1129  {
1130  i = first;
1131  j = last;
1132  kj = j-kind+1;
1133  while (j-i > tr)
1134  {
1135  if (i < cind)
1136  {
1137  alf = (ub-newkv[i])/(ua-newkv[i]);
1138  for (l = 0; l < size; l++)
1139  {
1140  newp.slice(i,l) = alf*newp.slice(i,l)-(1.0-alf)*newp.slice(i-1,l);
1141  }
1142  }
1143  if (j >= lbz)
1144  {
1145  if ((j-tr) <= (kind-ph+oldr))
1146  {
1147  gam = (ub-newkv[j-tr])/den;
1148  for (l = 0; l < size; l++)
1149  {
1150  ebpts(kj,l) = gam*ebpts(kj,l) + (1.0-gam)*ebpts(kj+1,l);
1151  }
1152  }
1153  else
1154  {
1155  for (l = 0; l < size; l++)
1156  {
1157  ebpts(kj,l) = bet*ebpts(kj,l) + (1.0-bet)*ebpts(kj+1,l);
1158  }
1159  }
1160  }
1161  i = i+1;
1162  j = j-1;
1163  kj = kj-1;
1164  }
1165  first--;
1166  last++;
1167  }
1168  }
1169
1170  if (a != p)
1171  {
1172  for (i = 0; i < (ph-oldr); i++)
1173  {
1174  newkv[kind] = ua;
1175  kind = kind+1;
1176  }
1177  }
1178  for (j = lbz; j <= rbz; j++)
1179  {
1180  for (l = 0; l < size; l++)
1181  {
1182  newp.slice(cind,l) = ebpts(j,l);
1183  }
1184  cind = cind +1;
1185  }
1186
1187  if (b < m)
1188  {
1189  for (j = 0; j <r; j++)
1190  for (l = 0; l < size; l++)
1191  {
1192  bpts(j,l) = nextbpts(j,l);
1193  }
1194
1195  for (j = r; j <= p; j++)
1196  for (l = 0; l < size; l++)
1197  {
1198  bpts(j,l) = oldp.slice(b-p+j,l);
1199  }
1200
1201  a = b;
1202  b = b+1;
1203  ua = ub;
1204  }
1205  else
1206  {
1207  for (i = 0; i <= ph; i++)
1208  {
1209  newkv[kind+i] = ub;
1210  }
1211  }
1212  }
1213  newkv.GetElements();
1214
1215  swap(newpatch);
1216 }
1217
1218 void NURBSPatch::FlipDirection(int dir)
1219 {
1220  int size = SetLoopDirection(dir);
1221
1222  for (int id = 0; id < nd/2; id++)
1223  for (int i = 0; i < size; i++)
1224  {
1225  Swap<double>((*this).slice(id,i), (*this).slice(nd-1-id,i));
1226  }
1227  kv[dir]->Flip();
1228 }
1229
1230 void NURBSPatch::SwapDirections(int dir1, int dir2)
1231 {
1232  if (abs(dir1-dir2) == 2)
1233  {
1234  mfem_error("NURBSPatch::SwapDirections :"
1235  " directions 0 and 2 are not supported!");
1236  }
1237
1238  Array<const KnotVector *> nkv(kv);
1239
1240  Swap<const KnotVector *>(nkv[dir1], nkv[dir2]);
1241  NURBSPatch *newpatch = new NURBSPatch(nkv, Dim);
1242
1243  int size = SetLoopDirection(dir1);
1244  newpatch->SetLoopDirection(dir2);
1245
1246  for (int id = 0; id < nd; id++)
1247  for (int i = 0; i < size; i++)
1248  {
1249  (*newpatch).slice(id,i) = (*this).slice(id,i);
1250  }
1251
1252  swap(newpatch);
1253 }
1254
1255 void NURBSPatch::Rotate(double angle, double n[])
1256 {
1257  if (Dim == 3)
1258  {
1259  Rotate2D(angle);
1260  }
1261  else
1262  {
1263  if (n == NULL)
1264  {
1265  mfem_error("NURBSPatch::Rotate : Specify an angle for a 3D rotation.");
1266  }
1267
1268  Rotate3D(n, angle);
1269  }
1270 }
1271
1272 void NURBSPatch::Get2DRotationMatrix(double angle, DenseMatrix &T)
1273 {
1274  double s = sin(angle);
1275  double c = cos(angle);
1276
1277  T.SetSize(2);
1278  T(0,0) = c;
1279  T(0,1) = -s;
1280  T(1,0) = s;
1281  T(1,1) = c;
1282 }
1283
1284 void NURBSPatch::Rotate2D(double angle)
1285 {
1286  if (Dim != 3)
1287  {
1288  mfem_error("NURBSPatch::Rotate2D : not a NURBSPatch in 2D!");
1289  }
1290
1291  DenseMatrix T(2);
1292  Vector x(2), y(NULL, 2);
1293
1294  Get2DRotationMatrix(angle, T);
1295
1296  int size = 1;
1297  for (int i = 0; i < kv.Size(); i++)
1298  {
1299  size *= kv[i]->GetNCP();
1300  }
1301
1302  for (int i = 0; i < size; i++)
1303  {
1304  y.SetData(data + i*Dim);
1305  x = y;
1306  T.Mult(x, y);
1307  }
1308 }
1309
1310 void NURBSPatch::Get3DRotationMatrix(double n[], double angle, double r,
1311  DenseMatrix &T)
1312 {
1313  double c, s, c1;
1314  double l2 = n[0]*n[0] + n[1]*n[1] + n[2]*n[2];
1315  double l = sqrt(l2);
1316
1317  if (fabs(angle) == M_PI_2)
1318  {
1319  s = r*copysign(1., angle);
1320  c = 0.;
1321  c1 = -1.;
1322  }
1323  else if (fabs(angle) == M_PI)
1324  {
1325  s = 0.;
1326  c = -r;
1327  c1 = c - 1.;
1328  }
1329  else
1330  {
1331  s = r*sin(angle);
1332  c = r*cos(angle);
1333  c1 = c - 1.;
1334  }
1335
1336  T.SetSize(3);
1337
1338  T(0,0) = (n[0]*n[0] + (n[1]*n[1] + n[2]*n[2])*c)/l2;
1339  T(0,1) = -(n[0]*n[1]*c1)/l2 - (n[2]*s)/l;
1340  T(0,2) = -(n[0]*n[2]*c1)/l2 + (n[1]*s)/l;
1341  T(1,0) = -(n[0]*n[1]*c1)/l2 + (n[2]*s)/l;
1342  T(1,1) = (n[1]*n[1] + (n[0]*n[0] + n[2]*n[2])*c)/l2;
1343  T(1,2) = -(n[1]*n[2]*c1)/l2 - (n[0]*s)/l;
1344  T(2,0) = -(n[0]*n[2]*c1)/l2 - (n[1]*s)/l;
1345  T(2,1) = -(n[1]*n[2]*c1)/l2 + (n[0]*s)/l;
1346  T(2,2) = (n[2]*n[2] + (n[0]*n[0] + n[1]*n[1])*c)/l2;
1347 }
1348
1349 void NURBSPatch::Rotate3D(double n[], double angle)
1350 {
1351  if (Dim != 4)
1352  {
1353  mfem_error("NURBSPatch::Rotate3D : not a NURBSPatch in 3D!");
1354  }
1355
1356  DenseMatrix T(3);
1357  Vector x(3), y(NULL, 3);
1358
1359  Get3DRotationMatrix(n, angle, 1., T);
1360
1361  int size = 1;
1362  for (int i = 0; i < kv.Size(); i++)
1363  {
1364  size *= kv[i]->GetNCP();
1365  }
1366
1367  for (int i = 0; i < size; i++)
1368  {
1369  y.SetData(data + i*Dim);
1370  x = y;
1371  T.Mult(x, y);
1372  }
1373 }
1374
1375 int NURBSPatch::MakeUniformDegree(int degree)
1376 {
1377  int maxd = degree;
1378
1379  if (maxd == -1)
1380  {
1381  for (int dir = 0; dir < kv.Size(); dir++)
1382  {
1383  maxd = std::max(maxd, kv[dir]->GetOrder());
1384  }
1385  }
1386
1387  for (int dir = 0; dir < kv.Size(); dir++)
1388  {
1389  if (maxd > kv[dir]->GetOrder())
1390  {
1391  DegreeElevate(dir, maxd - kv[dir]->GetOrder());
1392  }
1393  }
1394
1395  return maxd;
1396 }
1397
1398 NURBSPatch *Interpolate(NURBSPatch &p1, NURBSPatch &p2)
1399 {
1400  if (p1.kv.Size() != p2.kv.Size() || p1.Dim != p2.Dim)
1401  {
1402  mfem_error("Interpolate(NURBSPatch &, NURBSPatch &)");
1403  }
1404
1405  int size = 1, dim = p1.Dim;
1406  Array<const KnotVector *> kv(p1.kv.Size() + 1);
1407
1408  for (int i = 0; i < p1.kv.Size(); i++)
1409  {
1410  if (p1.kv[i]->GetOrder() < p2.kv[i]->GetOrder())
1411  {
1412  p1.KnotInsert(i, *p2.kv[i]);
1413  p2.KnotInsert(i, *p1.kv[i]);
1414  }
1415  else
1416  {
1417  p2.KnotInsert(i, *p1.kv[i]);
1418  p1.KnotInsert(i, *p2.kv[i]);
1419  }
1420  kv[i] = p1.kv[i];
1421  size *= kv[i]->GetNCP();
1422  }
1423
1424  KnotVector &nkv = *(new KnotVector(1, 2));
1425  nkv[0] = nkv[1] = 0.0;
1426  nkv[2] = nkv[3] = 1.0;
1427  nkv.GetElements();
1428  kv.Last() = &nkv;
1429
1430  NURBSPatch *patch = new NURBSPatch(kv, dim);
1431  delete kv.Last();
1432
1433  for (int i = 0; i < size; i++)
1434  {
1435  for (int d = 0; d < dim; d++)
1436  {
1437  patch->data[i*dim+d] = p1.data[i*dim+d];
1438  patch->data[(i+size)*dim+d] = p2.data[i*dim+d];
1439  }
1440  }
1441
1442  return patch;
1443 }
1444
1445 NURBSPatch *Revolve3D(NURBSPatch &patch, double n[], double ang, int times)
1446 {
1447  if (patch.Dim != 4)
1448  {
1449  mfem_error("Revolve3D(NURBSPatch &, double [], double)");
1450  }
1451
1452  int size = 1, ns;
1453  Array<const KnotVector *> nkv(patch.kv.Size() + 1);
1454
1455  for (int i = 0; i < patch.kv.Size(); i++)
1456  {
1457  nkv[i] = patch.kv[i];
1458  size *= nkv[i]->GetNCP();
1459  }
1460  ns = 2*times + 1;
1461  KnotVector &lkv = *(new KnotVector(2, ns));
1462  nkv.Last() = &lkv;
1463  lkv[0] = lkv[1] = lkv[2] = 0.0;
1464  for (int i = 1; i < times; i++)
1465  {
1466  lkv[2*i+1] = lkv[2*i+2] = i;
1467  }
1468  lkv[ns] = lkv[ns+1] = lkv[ns+2] = times;
1469  lkv.GetElements();
1470  NURBSPatch *newpatch = new NURBSPatch(nkv, 4);
1471  delete nkv.Last();
1472
1473  DenseMatrix T(3), T2(3);
1474  Vector u(NULL, 3), v(NULL, 3);
1475
1476  NURBSPatch::Get3DRotationMatrix(n, ang, 1., T);
1477  double c = cos(ang/2);
1478  NURBSPatch::Get3DRotationMatrix(n, ang/2, 1./c, T2);
1479  T2 *= c;
1480
1481  double *op = patch.data, *np;
1482  for (int i = 0; i < size; i++)
1483  {
1484  np = newpatch->data + 4*i;
1485  for (int j = 0; j < 4; j++)
1486  {
1487  np[j] = op[j];
1488  }
1489  for (int j = 0; j < times; j++)
1490  {
1491  u.SetData(np);
1492  v.SetData(np += 4*size);
1493  T2.Mult(u, v);
1494  v[3] = c*u[3];
1495  v.SetData(np += 4*size);
1496  T.Mult(u, v);
1497  v[3] = u[3];
1498  }
1499  op += 4;
1500  }
1501
1502  return newpatch;
1503 }
1504
1505
1506 NURBSExtension::NURBSExtension(const NURBSExtension &orig)
1507  : mOrder(orig.mOrder), mOrders(orig.mOrders),
1508  NumOfKnotVectors(orig.NumOfKnotVectors),
1509  NumOfVertices(orig.NumOfVertices),
1510  NumOfElements(orig.NumOfElements),
1511  NumOfBdrElements(orig.NumOfBdrElements),
1512  NumOfDofs(orig.NumOfDofs),
1513  NumOfActiveVertices(orig.NumOfActiveVertices),
1514  NumOfActiveElems(orig.NumOfActiveElems),
1515  NumOfActiveBdrElems(orig.NumOfActiveBdrElems),
1516  NumOfActiveDofs(orig.NumOfActiveDofs),
1517  activeVert(orig.activeVert),
1518  activeElem(orig.activeElem),
1519  activeBdrElem(orig.activeBdrElem),
1520  activeDof(orig.activeDof),
1521  patchTopo(new Mesh(*orig.patchTopo)),
1522  own_topo(true),
1523  edge_to_knot(orig.edge_to_knot),
1524  knotVectors(orig.knotVectors.Size()), // knotVectors are copied in the body
1525  knotVectorsCompr(orig.knotVectorsCompr.Size()),
1526  weights(orig.weights),
1527  d_to_d(orig.d_to_d),
1528  master(orig.master),
1529  slave(orig.slave),
1530  v_meshOffsets(orig.v_meshOffsets),
1531  e_meshOffsets(orig.e_meshOffsets),
1532  f_meshOffsets(orig.f_meshOffsets),
1533  p_meshOffsets(orig.p_meshOffsets),
1534  v_spaceOffsets(orig.v_spaceOffsets),
1535  e_spaceOffsets(orig.e_spaceOffsets),
1536  f_spaceOffsets(orig.f_spaceOffsets),
1537  p_spaceOffsets(orig.p_spaceOffsets),
1538  el_dof(orig.el_dof ? new Table(*orig.el_dof) : NULL),
1539  bel_dof(orig.bel_dof ? new Table(*orig.bel_dof) : NULL),
1540  el_to_patch(orig.el_to_patch),
1541  bel_to_patch(orig.bel_to_patch),
1542  el_to_IJK(orig.el_to_IJK),
1543  bel_to_IJK(orig.bel_to_IJK),
1544  patches(orig.patches.Size()) // patches are copied in the body
1545 {
1546  // Copy the knot vectors:
1547  for (int i = 0; i < knotVectors.Size(); i++)
1548  {
1549  knotVectors[i] = new KnotVector(*orig.knotVectors[i]);
1550  }
1551  CreateComprehensiveKV();
1552
1553  // Copy the patches:
1554  for (int p = 0; p < patches.Size(); p++)
1555  {
1556  patches[p] = new NURBSPatch(*orig.patches[p]);
1557  }
1558 }
1559
1560 NURBSExtension::NURBSExtension(std::istream &input)
1561 {
1563  patchTopo = new Mesh;
1565  own_topo = 1;
1566
1567  CheckPatches();
1568  // CheckBdrPatches();
1569
1570  skip_comment_lines(input, '#');
1571
1572  // Read knotvectors or patches
1573  string ident;
1574  input >> ws >> ident; // 'knotvectors' or 'patches'
1575  if (ident == "knotvectors")
1576  {
1577  input >> NumOfKnotVectors;
1578  knotVectors.SetSize(NumOfKnotVectors);
1579  for (int i = 0; i < NumOfKnotVectors; i++)
1580  {
1581  knotVectors[i] = new KnotVector(input);
1582  }
1583  }
1584  else if (ident == "patches")
1585  {
1586  patches.SetSize(GetNP());
1587  for (int p = 0; p < patches.Size(); p++)
1588  {
1589  skip_comment_lines(input, '#');
1590  patches[p] = new NURBSPatch(input);
1591  }
1592
1593  NumOfKnotVectors = 0;
1594  for (int i = 0; i < patchTopo->GetNEdges(); i++)
1595  if (NumOfKnotVectors < KnotInd(i))
1596  {
1597  NumOfKnotVectors = KnotInd(i);
1598  }
1599  NumOfKnotVectors++;
1600  knotVectors.SetSize(NumOfKnotVectors);
1601  knotVectors = NULL;
1602
1603  Array<int> edges, oedge;
1604  for (int p = 0; p < patches.Size(); p++)
1605  {
1606  if (Dimension() == 1)
1607  {
1608  if (knotVectors[KnotInd(p)] == NULL)
1609  {
1610  knotVectors[KnotInd(p)] =
1611  new KnotVector(*patches[p]->GetKV(0));
1612  }
1613  }
1614  if (Dimension() == 2)
1615  {
1616  patchTopo->GetElementEdges(p, edges, oedge);
1617  if (knotVectors[KnotInd(edges[0])] == NULL)
1618  {
1619  knotVectors[KnotInd(edges[0])] =
1620  new KnotVector(*patches[p]->GetKV(0));
1621  }
1622  if (knotVectors[KnotInd(edges[1])] == NULL)
1623  {
1624  knotVectors[KnotInd(edges[1])] =
1625  new KnotVector(*patches[p]->GetKV(1));
1626  }
1627  }
1628  else if (Dimension() == 3)
1629  {
1630  patchTopo->GetElementEdges(p, edges, oedge);
1631  if (knotVectors[KnotInd(edges[0])] == NULL)
1632  {
1633  knotVectors[KnotInd(edges[0])] =
1634  new KnotVector(*patches[p]->GetKV(0));
1635  }
1636  if (knotVectors[KnotInd(edges[3])] == NULL)
1637  {
1638  knotVectors[KnotInd(edges[3])] =
1639  new KnotVector(*patches[p]->GetKV(1));
1640  }
1641  if (knotVectors[KnotInd(edges[8])] == NULL)
1642  {
1643  knotVectors[KnotInd(edges[8])] =
1644  new KnotVector(*patches[p]->GetKV(2));
1645  }
1646  }
1647  }
1648  }
1649  else
1650  {
1651  MFEM_ABORT("invalid section: " << ident);
1652  }
1653
1654  CreateComprehensiveKV();
1655
1656  SetOrdersFromKnotVectors();
1657
1658  GenerateOffsets();
1659  CountElements();
1660  CountBdrElements();
1661  // NumOfVertices, NumOfElements, NumOfBdrElements, NumOfDofs
1662
1663  skip_comment_lines(input, '#');
1664
1665  // Check for a list of mesh elements
1666  if (patches.Size() == 0)
1667  {
1668  input >> ws >> ident;
1669  }
1670  if (patches.Size() == 0 && ident == "mesh_elements")
1671  {
1672  input >> NumOfActiveElems;
1673  activeElem.SetSize(GetGNE());
1674  activeElem = false;
1675  int glob_elem;
1676  for (int i = 0; i < NumOfActiveElems; i++)
1677  {
1678  input >> glob_elem;
1679  activeElem[glob_elem] = true;
1680  }
1681
1682  skip_comment_lines(input, '#');
1683  input >> ws >> ident;
1684  }
1685  else
1686  {
1687  NumOfActiveElems = NumOfElements;
1688  activeElem.SetSize(NumOfElements);
1689  activeElem = true;
1690  }
1691
1692  GenerateActiveVertices();
1693  InitDofMap();
1694  GenerateElementDofTable();
1695  GenerateActiveBdrElems();
1696  GenerateBdrElementDofTable();
1697
1698  // periodic
1699  if (ident == "periodic")
1700  {
1703
1704  skip_comment_lines(input, '#');
1705  input >> ws >> ident;
1706  }
1707
1708  if (patches.Size() == 0)
1709  {
1710  // weights
1711  if (ident == "weights")
1712  {
1714  }
1715  else // e.g. ident = "unitweights" or "autoweights"
1716  {
1717  weights.SetSize(GetNDof());
1718  weights = 1.0;
1719  }
1720  }
1721
1722  // periodic
1723  ConnectBoundaries();
1724 }
1725
1726 NURBSExtension::NURBSExtension(NURBSExtension *parent, int newOrder)
1727 {
1728  patchTopo = parent->patchTopo;
1729  own_topo = 0;
1730
1731  parent->edge_to_knot.Copy(edge_to_knot);
1732
1733  NumOfKnotVectors = parent->GetNKV();
1734  knotVectors.SetSize(NumOfKnotVectors);
1735  knotVectorsCompr.SetSize(parent->GetNP()*parent->Dimension());
1736  const Array<int> &pOrders = parent->GetOrders();
1737  for (int i = 0; i < NumOfKnotVectors; i++)
1738  {
1739  if (newOrder > pOrders[i])
1740  {
1741  knotVectors[i] =
1742  parent->GetKnotVector(i)->DegreeElevate(newOrder - pOrders[i]);
1743  }
1744  else
1745  {
1746  knotVectors[i] = new KnotVector(*parent->GetKnotVector(i));
1747  }
1748  }
1749  CreateComprehensiveKV();
1750
1751  // copy some data from parent
1752  NumOfElements = parent->NumOfElements;
1753  NumOfBdrElements = parent->NumOfBdrElements;
1754
1755  SetOrdersFromKnotVectors();
1756
1757  GenerateOffsets(); // dof offsets will be different from parent
1758
1759  NumOfActiveVertices = parent->NumOfActiveVertices;
1760  NumOfActiveElems = parent->NumOfActiveElems;
1761  NumOfActiveBdrElems = parent->NumOfActiveBdrElems;
1762  parent->activeVert.Copy(activeVert);
1763  InitDofMap();
1764  parent->activeElem.Copy(activeElem);
1765  parent->activeBdrElem.Copy(activeBdrElem);
1766
1767  GenerateElementDofTable();
1768  GenerateBdrElementDofTable();
1769
1770  weights.SetSize(GetNDof());
1771  weights = 1.0;
1772
1773  // periodic
1774  parent->master.Copy(master);
1775  parent->slave.Copy(slave);
1776  ConnectBoundaries();
1777 }
1778
1779 NURBSExtension::NURBSExtension(NURBSExtension *parent,
1780  const Array<int> &newOrders)
1781 {
1782  newOrders.Copy(mOrders);
1783  SetOrderFromOrders();
1784
1785  patchTopo = parent->patchTopo;
1786  own_topo = 0;
1787
1788  parent->edge_to_knot.Copy(edge_to_knot);
1789
1790  NumOfKnotVectors = parent->GetNKV();
1791  MFEM_VERIFY(mOrders.Size() == NumOfKnotVectors, "invalid newOrders array");
1792  knotVectors.SetSize(NumOfKnotVectors);
1793  const Array<int> &pOrders = parent->GetOrders();
1794
1795  for (int i = 0; i < NumOfKnotVectors; i++)
1796  {
1797  if (mOrders[i] > pOrders[i])
1798  {
1799  knotVectors[i] =
1800  parent->GetKnotVector(i)->DegreeElevate(mOrders[i] - pOrders[i]);
1801  }
1802  else
1803  {
1804  knotVectors[i] = new KnotVector(*parent->GetKnotVector(i));
1805  }
1806  }
1807  CreateComprehensiveKV();
1808
1809  // copy some data from parent
1810  NumOfElements = parent->NumOfElements;
1811  NumOfBdrElements = parent->NumOfBdrElements;
1812
1813  GenerateOffsets(); // dof offsets will be different from parent
1814
1815  NumOfActiveVertices = parent->NumOfActiveVertices;
1816  NumOfActiveElems = parent->NumOfActiveElems;
1817  NumOfActiveBdrElems = parent->NumOfActiveBdrElems;
1818  parent->activeVert.Copy(activeVert);
1819  InitDofMap();
1820  parent->activeElem.Copy(activeElem);
1821  parent->activeBdrElem.Copy(activeBdrElem);
1822
1823  GenerateElementDofTable();
1824  GenerateBdrElementDofTable();
1825
1826  weights.SetSize(GetNDof());
1827  weights = 1.0;
1828
1829  parent->master.Copy(master);
1830  parent->slave.Copy(slave);
1831  ConnectBoundaries();
1832 }
1833
1834 NURBSExtension::NURBSExtension(Mesh *mesh_array[], int num_pieces)
1835 {
1836  NURBSExtension *parent = mesh_array[0]->NURBSext;
1837
1838  if (!parent->own_topo)
1839  {
1840  mfem_error("NURBSExtension::NURBSExtension :\n"
1841  " parent does not own the patch topology!");
1842  }
1843  patchTopo = parent->patchTopo;
1844  own_topo = 1;
1845  parent->own_topo = 0;
1846
1847  parent->edge_to_knot.Copy(edge_to_knot);
1848
1849  parent->GetOrders().Copy(mOrders);
1850  mOrder = parent->GetOrder();
1851
1852  NumOfKnotVectors = parent->GetNKV();
1853  knotVectors.SetSize(NumOfKnotVectors);
1854  for (int i = 0; i < NumOfKnotVectors; i++)
1855  {
1856  knotVectors[i] = new KnotVector(*parent->GetKnotVector(i));
1857  }
1858  CreateComprehensiveKV();
1859
1860  GenerateOffsets();
1861  CountElements();
1862  CountBdrElements();
1863
1864  // assuming the meshes define a partitioning of all the elements
1865  NumOfActiveElems = NumOfElements;
1866  activeElem.SetSize(NumOfElements);
1867  activeElem = true;
1868
1869  GenerateActiveVertices();
1870  InitDofMap();
1871  GenerateElementDofTable();
1872  GenerateActiveBdrElems();
1873  GenerateBdrElementDofTable();
1874
1875  weights.SetSize(GetNDof());
1876  MergeWeights(mesh_array, num_pieces);
1877 }
1878
1879 NURBSExtension::~NURBSExtension()
1880 {
1881  if (patches.Size() == 0)
1882  {
1883  delete bel_dof;
1884  delete el_dof;
1885  }
1886
1887  for (int i = 0; i < knotVectors.Size(); i++)
1888  {
1889  delete knotVectors[i];
1890  }
1891
1892  for (int i = 0; i < knotVectorsCompr.Size(); i++)
1893  {
1894  delete knotVectorsCompr[i];
1895  }
1896
1897  for (int i = 0; i < patches.Size(); i++)
1898  {
1899  delete patches[i];
1900  }
1901
1902  if (own_topo)
1903  {
1904  delete patchTopo;
1905  }
1906 }
1907
1908 void NURBSExtension::Print(std::ostream &os) const
1909 {
1910  patchTopo->PrintTopo(os, edge_to_knot);
1911  if (patches.Size() == 0)
1912  {
1913  os << "\nknotvectors\n" << NumOfKnotVectors << '\n';
1914  for (int i = 0; i < NumOfKnotVectors; i++)
1915  {
1916  knotVectors[i]->Print(os);
1917  }
1918
1919  if (NumOfActiveElems < NumOfElements)
1920  {
1921  os << "\nmesh_elements\n" << NumOfActiveElems << '\n';
1922  for (int i = 0; i < NumOfElements; i++)
1923  if (activeElem[i])
1924  {
1925  os << i << '\n';
1926  }
1927  }
1928
1929  os << "\nweights\n";
1930  weights.Print(os, 1);
1931  }
1932  else
1933  {
1934  os << "\npatches\n";
1935  for (int p = 0; p < patches.Size(); p++)
1936  {
1937  os << "\n# patch " << p << "\n\n";
1938  patches[p]->Print(os);
1939  }
1940  }
1941 }
1942
1943 void NURBSExtension::PrintCharacteristics(std::ostream &os) const
1944 {
1945  os <<
1946  "NURBS Mesh entity sizes:\n"
1947  "Dimension = " << Dimension() << "\n"
1948  "Unique Orders = ";
1949  Array<int> unique_orders(mOrders);
1950  unique_orders.Sort();
1951  unique_orders.Unique();
1952  unique_orders.Print(os, unique_orders.Size());
1953  os <<
1954  "NumOfKnotVectors = " << GetNKV() << "\n"
1955  "NumOfPatches = " << GetNP() << "\n"
1956  "NumOfBdrPatches = " << GetNBP() << "\n"
1957  "NumOfVertices = " << GetGNV() << "\n"
1958  "NumOfElements = " << GetGNE() << "\n"
1959  "NumOfBdrElements = " << GetGNBE() << "\n"
1960  "NumOfDofs = " << GetNTotalDof() << "\n"
1961  "NumOfActiveVertices = " << GetNV() << "\n"
1962  "NumOfActiveElems = " << GetNE() << "\n"
1963  "NumOfActiveBdrElems = " << GetNBE() << "\n"
1964  "NumOfActiveDofs = " << GetNDof() << '\n';
1965  for (int i = 0; i < NumOfKnotVectors; i++)
1966  {
1967  os << ' ' << i + 1 << ") ";
1968  knotVectors[i]->Print(os);
1969  }
1970  os << endl;
1971 }
1972
1973 void NURBSExtension::PrintFunctions(const char *basename, int samples) const
1974 {
1975  std::ofstream os;
1976  for (int i = 0; i < NumOfKnotVectors; i++)
1977  {
1978  std::ostringstream filename;
1979  filename << basename<<"_"<<i<<".dat";
1980  os.open(filename.str().c_str());
1981  knotVectors[i]->PrintFunctions(os,samples);
1982  os.close();
1983  }
1984 }
1985
1986 void NURBSExtension::InitDofMap()
1987 {
1988  master.SetSize(0);
1989  slave.SetSize(0);
1990  d_to_d.SetSize(0);
1991 }
1992
1993 void NURBSExtension::ConnectBoundaries(Array<int> &bnds0, Array<int> &bnds1)
1994 {
1995  bnds0.Copy(master);
1996  bnds1.Copy(slave);
1997  ConnectBoundaries();
1998 }
1999
2000 void NURBSExtension::ConnectBoundaries()
2001 {
2002  if (master.Size() != slave.Size())
2003  {
2004  mfem_error("NURBSExtension::ConnectBoundaries() boundary lists not of equal size");
2005  }
2006  if (master.Size() == 0 ) { return; }
2007
2008  // Initialize d_to_d
2009  d_to_d.SetSize(NumOfDofs);
2010  for (int i = 0; i < NumOfDofs; i++) { d_to_d[i] = i; }
2011
2012  // Connect
2013  for (int i = 0; i < master.Size(); i++)
2014  {
2015  int bnd0 = -1, bnd1 = -1;
2016  for (int b = 0; b < GetNBP(); b++)
2017  {
2018  if (master[i] == patchTopo->GetBdrAttribute(b)) { bnd0 = b; }
2019  if (slave[i]== patchTopo->GetBdrAttribute(b)) { bnd1 = b; }
2020  }
2023
2024  if (Dimension() == 1)
2025  {
2026  ConnectBoundaries1D(bnd0, bnd1);
2027  }
2028  else if (Dimension() == 2)
2029  {
2030  ConnectBoundaries2D(bnd0, bnd1);
2031  }
2032  else
2033  {
2034  ConnectBoundaries3D(bnd0, bnd1);
2035  }
2036  }
2037
2038  // Clean d_to_d
2039  Array<int> tmp(d_to_d.Size()+1);
2040  tmp = 0;
2041
2042  for (int i = 0; i < d_to_d.Size(); i++)
2043  {
2044  tmp[d_to_d[i]] = 1;
2045  }
2046
2047  int cnt = 0;
2048  for (int i = 0; i < tmp.Size(); i++)
2049  {
2050  if (tmp[i] == 1) { tmp[i] = cnt++; }
2051  }
2052  NumOfDofs = cnt;
2053
2054  for (int i = 0; i < d_to_d.Size(); i++)
2055  {
2056  d_to_d[i] = tmp[d_to_d[i]];
2057  }
2058
2059  // Finalize
2060  if (el_dof) { delete el_dof; }
2061  if (bel_dof) { delete bel_dof; }
2062  GenerateElementDofTable();
2063  GenerateBdrElementDofTable();
2064 }
2065
2066 void NURBSExtension::ConnectBoundaries1D(int bnd0, int bnd1)
2067 {
2068  NURBSPatchMap p2g0(this);
2069  NURBSPatchMap p2g1(this);
2070
2071  int okv0[1],okv1[1];
2072  const KnotVector *kv0[1],*kv1[1];
2073
2074  p2g0.SetBdrPatchDofMap(bnd0, kv0, okv0);
2075  p2g1.SetBdrPatchDofMap(bnd1, kv1, okv1);
2076
2077  d_to_d[p2g0(0)] = d_to_d[p2g1(0)];
2078 }
2079
2080 void NURBSExtension::ConnectBoundaries2D(int bnd0, int bnd1)
2081 {
2082  NURBSPatchMap p2g0(this);
2083  NURBSPatchMap p2g1(this);
2084
2085  int okv0[1],okv1[1];
2086  const KnotVector *kv0[1],*kv1[1];
2087
2088  p2g0.SetBdrPatchDofMap(bnd0, kv0, okv0);
2089  p2g1.SetBdrPatchDofMap(bnd1, kv1, okv1);
2090
2091  int nx = p2g0.nx();
2092  int nks0 = kv0[0]->GetNKS();
2093
2094 #ifdef MFEM_DEBUG
2095  bool compatible = true;
2096  if (p2g0.nx() != p2g1.nx()) { compatible = false; }
2097  if (kv0[0]->GetNKS() != kv1[0]->GetNKS()) { compatible = false; }
2098  if (kv0[0]->GetOrder() != kv1[0]->GetOrder()) { compatible = false; }
2099
2100  if (!compatible)
2101  {
2102  mfem::out<<p2g0.nx()<<" "<<p2g1.nx()<<endl;
2103  mfem::out<<kv0[0]->GetNKS()<<" "<<kv1[0]->GetNKS()<<endl;
2104  mfem::out<<kv0[0]->GetOrder()<<" "<<kv1[0]->GetOrder()<<endl;
2105  mfem_error("NURBS boundaries not compatible");
2106  }
2107 #endif
2108
2109  for (int i = 0; i < nks0; i++)
2110  {
2111  if (kv0[0]->isElement(i))
2112  {
2113  if (!kv1[0]->isElement(i)) { mfem_error("isElement does not match"); }
2114  for (int ii = 0; ii <= kv0[0]->GetOrder(); ii++)
2115  {
2116  int ii0 = (okv0[0] >= 0) ? (i+ii) : (nx-i-ii);
2117  int ii1 = (okv1[0] >= 0) ? (i+ii) : (nx-i-ii);
2118
2119  d_to_d[p2g0(ii0)] = d_to_d[p2g1(ii1)];
2120  }
2121
2122  }
2123  }
2124 }
2125
2126 void NURBSExtension::ConnectBoundaries3D(int bnd0, int bnd1)
2127 {
2128  NURBSPatchMap p2g0(this);
2129  NURBSPatchMap p2g1(this);
2130
2131  int okv0[2],okv1[2];
2132  const KnotVector *kv0[2],*kv1[2];
2133
2134  p2g0.SetBdrPatchDofMap(bnd0, kv0, okv0);
2135  p2g1.SetBdrPatchDofMap(bnd1, kv1, okv1);
2136
2137  int nx = p2g0.nx();
2138  int ny = p2g0.ny();
2139
2140  int nks0 = kv0[0]->GetNKS();
2141  int nks1 = kv0[1]->GetNKS();
2142
2143 #ifdef MFEM_DEBUG
2144  bool compatible = true;
2145  if (p2g0.nx() != p2g1.nx()) { compatible = false; }
2146  if (p2g0.ny() != p2g1.ny()) { compatible = false; }
2147
2148  if (kv0[0]->GetNKS() != kv1[0]->GetNKS()) { compatible = false; }
2149  if (kv0[1]->GetNKS() != kv1[1]->GetNKS()) { compatible = false; }
2150
2151  if (kv0[0]->GetOrder() != kv1[0]->GetOrder()) { compatible = false; }
2152  if (kv0[1]->GetOrder() != kv1[1]->GetOrder()) { compatible = false; }
2153
2154  if (!compatible)
2155  {
2156  mfem::out<<p2g0.nx()<<" "<<p2g1.nx()<<endl;
2157  mfem::out<<p2g0.ny()<<" "<<p2g1.ny()<<endl;
2158
2159  mfem::out<<kv0[0]->GetNKS()<<" "<<kv1[0]->GetNKS()<<endl;
2160  mfem::out<<kv0[1]->GetNKS()<<" "<<kv1[1]->GetNKS()<<endl;
2161
2162  mfem::out<<kv0[0]->GetOrder()<<" "<<kv1[0]->GetOrder()<<endl;
2163  mfem::out<<kv0[1]->GetOrder()<<" "<<kv1[1]->GetOrder()<<endl;
2164  mfem_error("NURBS boundaries not compatible");
2165  }
2166 #endif
2167
2168  for (int j = 0; j < nks1; j++)
2169  {
2170  if (kv0[1]->isElement(j))
2171  {
2172  if (!kv1[1]->isElement(j)) { mfem_error("isElement does not match #1"); }
2173  for (int i = 0; i < nks0; i++)
2174  {
2175  if (kv0[0]->isElement(i))
2176  {
2177  if (!kv1[0]->isElement(i)) { mfem_error("isElement does not match #0"); }
2178  for (int jj = 0; jj <= kv0[1]->GetOrder(); jj++)
2179  {
2180  int jj0 = (okv0[1] >= 0) ? (j+jj) : (ny-j-jj);
2181  int jj1 = (okv1[1] >= 0) ? (j+jj) : (ny-j-jj);
2182
2183  for (int ii = 0; ii <= kv0[0]->GetOrder(); ii++)
2184  {
2185  int ii0 = (okv0[0] >= 0) ? (i+ii) : (nx-i-ii);
2186  int ii1 = (okv1[0] >= 0) ? (i+ii) : (nx-i-ii);
2187
2188  d_to_d[p2g0(ii0,jj0)] = d_to_d[p2g1(ii1,jj1)];
2189  }
2190  }
2191  }
2192  }
2193  }
2194  }
2195 }
2196
2198 {
2199  int vert[8], nv, g_el, nx, ny, nz, dim = Dimension();
2200
2201  NURBSPatchMap p2g(this);
2202  const KnotVector *kv[3];
2203
2204  g_el = 0;
2205  activeVert.SetSize(GetGNV());
2206  activeVert = -1;
2207  for (int p = 0; p < GetNP(); p++)
2208  {
2209  p2g.SetPatchVertexMap(p, kv);
2210
2211  nx = p2g.nx();
2212  ny = (dim >= 2) ? p2g.ny() : 1;
2213  nz = (dim == 3) ? p2g.nz() : 1;
2214
2215  for (int k = 0; k < nz; k++)
2216  {
2217  for (int j = 0; j < ny; j++)
2218  {
2219  for (int i = 0; i < nx; i++)
2220  {
2221  if (activeElem[g_el])
2222  {
2223  if (dim == 1)
2224  {
2225  vert[0] = p2g(i );
2226  vert[1] = p2g(i+1);
2227  nv = 2;
2228  }
2229  else if (dim == 2)
2230  {
2231  vert[0] = p2g(i, j );
2232  vert[1] = p2g(i+1,j );
2233  vert[2] = p2g(i+1,j+1);
2234  vert[3] = p2g(i, j+1);
2235  nv = 4;
2236  }
2237  else
2238  {
2239  vert[0] = p2g(i, j, k);
2240  vert[1] = p2g(i+1,j, k);
2241  vert[2] = p2g(i+1,j+1,k);
2242  vert[3] = p2g(i, j+1,k);
2243
2244  vert[4] = p2g(i, j, k+1);
2245  vert[5] = p2g(i+1,j, k+1);
2246  vert[6] = p2g(i+1,j+1,k+1);
2247  vert[7] = p2g(i, j+1,k+1);
2248  nv = 8;
2249  }
2250
2251  for (int v = 0; v < nv; v++)
2252  {
2253  activeVert[vert[v]] = 1;
2254  }
2255  }
2256  g_el++;
2257  }
2258  }
2259  }
2260  }
2261
2262  NumOfActiveVertices = 0;
2263  for (int i = 0; i < GetGNV(); i++)
2264  if (activeVert[i] == 1)
2265  {
2266  activeVert[i] = NumOfActiveVertices++;
2267  }
2268 }
2269
2271 {
2272  int dim = Dimension();
2274
2275  activeBdrElem.SetSize(GetGNBE());
2276  if (GetGNE() == GetNE())
2277  {
2278  activeBdrElem = true;
2279  NumOfActiveBdrElems = GetGNBE();
2280  return;
2281  }
2282  activeBdrElem = false;
2283  NumOfActiveBdrElems = 0;
2284  // the mesh will generate the actual boundary including boundary
2285  // elements that are not on boundary patches. we use this for
2286  // visualization of processor boundaries
2287
2288  // TODO: generate actual boundary?
2289 }
2290
2291
2292 void NURBSExtension::MergeWeights(Mesh *mesh_array[], int num_pieces)
2293 {
2294  Array<int> lelem_elem;
2295
2296  for (int i = 0; i < num_pieces; i++)
2297  {
2298  NURBSExtension *lext = mesh_array[i]->NURBSext;
2299
2300  lext->GetElementLocalToGlobal(lelem_elem);
2301
2302  for (int lel = 0; lel < lext->GetNE(); lel++)
2303  {
2304  int gel = lelem_elem[lel];
2305
2306  int nd = el_dof->RowSize(gel);
2307  int *gdofs = el_dof->GetRow(gel);
2308  int *ldofs = lext->el_dof->GetRow(lel);
2309  for (int j = 0; j < nd; j++)
2310  {
2311  weights(gdofs[j]) = lext->weights(ldofs[j]);
2312  }
2313  }
2314  }
2315 }
2316
2318  GridFunction *gf_array[], int num_pieces, GridFunction &merged)
2319 {
2320  FiniteElementSpace *gfes = merged.FESpace();
2321  Array<int> lelem_elem, dofs;
2322  Vector lvec;
2323
2324  for (int i = 0; i < num_pieces; i++)
2325  {
2326  FiniteElementSpace *lfes = gf_array[i]->FESpace();
2327  NURBSExtension *lext = lfes->GetMesh()->NURBSext;
2328
2329  lext->GetElementLocalToGlobal(lelem_elem);
2330
2331  for (int lel = 0; lel < lext->GetNE(); lel++)
2332  {
2333  lfes->GetElementVDofs(lel, dofs);
2334  gf_array[i]->GetSubVector(dofs, lvec);
2335
2336  gfes->GetElementVDofs(lelem_elem[lel], dofs);
2337  merged.SetSubVector(dofs, lvec);
2338  }
2339  }
2340 }
2341
2343 {
2344  if (Dimension() == 1 ) { return; }
2345
2346  Array<int> edges;
2347  Array<int> oedge;
2348
2349  for (int p = 0; p < GetNP(); p++)
2350  {
2351  patchTopo->GetElementEdges(p, edges, oedge);
2352
2353  for (int i = 0; i < edges.Size(); i++)
2354  {
2355  edges[i] = edge_to_knot[edges[i]];
2356  if (oedge[i] < 0)
2357  {
2358  edges[i] = -1 - edges[i];
2359  }
2360  }
2361
2362  if ((Dimension() == 2 &&
2363  (edges[0] != -1 - edges[2] || edges[1] != -1 - edges[3])) ||
2364
2365  (Dimension() == 3 &&
2366  (edges[0] != edges[2] || edges[0] != edges[4] ||
2367  edges[0] != edges[6] || edges[1] != edges[3] ||
2368  edges[1] != edges[5] || edges[1] != edges[7] ||
2369  edges[8] != edges[9] || edges[8] != edges[10] ||
2370  edges[8] != edges[11])))
2371  {
2372  mfem::err << "NURBSExtension::CheckPatch (patch = " << p
2373  << ")\n Inconsistent edge-to-knot mapping!\n";
2374  mfem_error();
2375  }
2376  }
2377 }
2378
2380 {
2381  Array<int> edges;
2382  Array<int> oedge;
2383
2384  for (int p = 0; p < GetNBP(); p++)
2385  {
2386  patchTopo->GetBdrElementEdges(p, edges, oedge);
2387
2388  for (int i = 0; i < edges.Size(); i++)
2389  {
2390  edges[i] = edge_to_knot[edges[i]];
2391  if (oedge[i] < 0)
2392  {
2393  edges[i] = -1 - edges[i];
2394  }
2395  }
2396
2397  if ((Dimension() == 2 && (edges[0] < 0)) ||
2398  (Dimension() == 3 && (edges[0] < 0 || edges[1] < 0)))
2399  {
2400  mfem::err << "NURBSExtension::CheckBdrPatch (boundary patch = "
2401  << p << ") : Bad orientation!\n";
2402  mfem_error();
2403  }
2404  }
2405 }
2406
2408 {
2409  // patchTopo->GetElementEdges is not yet implemented for 1D
2410  MFEM_VERIFY(Dimension()>1, "1D not yet implemented.");
2411
2412  kvdir.SetSize(Dimension());
2413  kvdir = 0;
2414
2415  Array<int> patchvert, edges, orient, edgevert;
2416
2417  patchTopo->GetElementVertices(p, patchvert);
2418
2419  patchTopo->GetElementEdges(p, edges, orient);
2420
2421  // Compare the vertices of the patches with the vertices of the knotvectors of knot2dge
2422  // Based on the match the orientation will be a 1 or a -1
2423  // -1: direction is flipped
2424  // 1: direction is not flipped
2425
2426
2427  for (int i = 0; i < edges.Size(); i++)
2428  {
2429  // First side
2430  patchTopo->GetEdgeVertices(edges[i], edgevert);
2431  if (edgevert[0] == patchvert[0] && edgevert[1] == patchvert[1])
2432  {
2433  kvdir[0] = 1;
2434  }
2435
2436  if (edgevert[0] == patchvert[1] && edgevert[1] == patchvert[0])
2437  {
2438  kvdir[0] = -1;
2439  }
2440
2441  // Second side
2442  if (edgevert[0] == patchvert[1] && edgevert[1] == patchvert[2])
2443  {
2444  kvdir[1] = 1;
2445  }
2446
2447  if (edgevert[0] == patchvert[2] && edgevert[1] == patchvert[1])
2448  {
2449  kvdir[1] = -1;
2450  }
2451  }
2452
2453  if (Dimension() == 3)
2454  {
2455  // Third side
2456  for (int i = 0; i < edges.Size(); i++)
2457  {
2458  patchTopo->GetEdgeVertices(edges[i], edgevert);
2459
2460  if (edgevert[0] == patchvert[0] && edgevert[1] == patchvert[4])
2461  {
2462  kvdir[2] = 1;
2463  }
2464
2465  if (edgevert[0] == patchvert[4] && edgevert[1] == patchvert[0])
2466  {
2467  kvdir[2] = -1;
2468  }
2469  }
2470  }
2471
2472  MFEM_VERIFY(kvdir.Find(0) == -1, "Could not find direction of knotvector.");
2473 }
2474
2476 {
2477  Array<int> edges, orient, kvdir;
2478  Array<int> e(Dimension());
2479
2480  // 1D: comprehensive and unique KV are the same
2481  if (Dimension() == 1)
2482  {
2483  knotVectorsCompr.SetSize(GetNKV());
2484  for (int i = 0; i < GetNKV(); i++)
2485  {
2486  knotVectorsCompr[i] = new KnotVector(*(KnotVec(i)));
2487  }
2488  return;
2489  }
2490  else if (Dimension() == 2)
2491  {
2492  knotVectorsCompr.SetSize(GetNP()*Dimension());
2493  e[0] = 0;
2494  e[1] = 1;
2495  }
2496  else if (Dimension() == 3)
2497  {
2498  knotVectorsCompr.SetSize(GetNP()*Dimension());
2499  e[0] = 0;
2500  e[1] = 3;
2501  e[2] = 8;
2502  }
2503
2504  for (int p = 0; p < GetNP(); p++)
2505  {
2506  CheckKVDirection(p, kvdir);
2507
2508  patchTopo->GetElementEdges(p, edges, orient);
2509
2510  for (int d = 0; d < Dimension(); d++)
2511  {
2512  // Indices in unique and comprehensive sets of the KnotVector
2513  int iun = edges[e[d]];
2514  int icomp = Dimension()*p+d;
2515
2516  knotVectorsCompr[icomp] = new KnotVector(*(KnotVec(iun)));
2517
2518  if (kvdir[d] == -1) {knotVectorsCompr[icomp]->Flip();}
2519  }
2520  }
2521
2522  MFEM_VERIFY(ConsistentKVSets(), "Mismatch in KnotVectors");
2523 }
2524
2526 {
2527  Array<int> e(Dimension());
2528
2529  // 1D: comprehensive and unique KV are the same
2530  if (Dimension() == 1)
2531  {
2532  for (int i = 0; i < GetNKV(); i++)
2533  {
2534  *(KnotVec(i)) = *(knotVectorsCompr[i]);
2535  }
2536  return;
2537  }
2538  else if (Dimension() == 2)
2539  {
2540  e[0] = 0;
2541  e[1] = 1;
2542  }
2543  else if (Dimension() == 3)
2544  {
2545  e[0] = 0;
2546  e[1] = 3;
2547  e[2] = 8;
2548  }
2549
2550  for (int p = 0; p < GetNP(); p++)
2551  {
2552  Array<int> edges, orient, kvdir;
2553
2554  patchTopo->GetElementEdges(p, edges, orient);
2555  CheckKVDirection(p, kvdir);
2556
2557  for ( int d = 0; d < Dimension(); d++)
2558  {
2559  bool flip = false;
2560  if (kvdir[d] == -1) {flip = true;}
2561
2562  // Indices in unique and comprehensive sets of the KnotVector
2563  int iun = edges[e[d]];
2564  int icomp = Dimension()*p+d;
2565
2566  // Check if difference in order
2567  int o1 = KnotVec(iun)->GetOrder();
2568  int o2 = knotVectorsCompr[icomp]->GetOrder();
2569  int diffo = abs(o1 - o2);
2570
2571  if (diffo)
2572  {
2573  // Update reduced set of knotvectors
2574  *(KnotVec(iun)) = *(knotVectorsCompr[icomp]);
2575
2576  // Give correct direction to unique knotvector.
2577  if (flip) { KnotVec(iun)->Flip(); }
2578  }
2579
2580  // Check if difference between knots
2581  Vector diffknot;
2582
2583  if (flip) { knotVectorsCompr[icomp]->Flip(); }
2584
2585  KnotVec(iun)->Difference(*(knotVectorsCompr[icomp]), diffknot);
2586
2587  if (flip) { knotVectorsCompr[icomp]->Flip(); }
2588
2589  if (diffknot.Size() > 0)
2590  {
2591  // Update reduced set of knotvectors
2592  *(KnotVec(iun)) = *(knotVectorsCompr[icomp]);
2593
2594  // Give correct direction to unique knotvector.
2595  if (flip) {KnotVec(iun)->Flip();}
2596  }
2597  }
2598  }
2599
2600  MFEM_VERIFY(ConsistentKVSets(), "Mismatch in KnotVectors");
2601 }
2602
2604 {
2605  // patchTopo->GetElementEdges is not yet implemented for 1D
2606  MFEM_VERIFY(Dimension()>1, "1D not yet implemented.");
2607
2608  Array<int> edges, orient, kvdir;
2609  Vector diff;
2610
2611  Array<int>e(Dimension());
2612
2613  e[0] = 0;
2614
2615  if (Dimension() == 2)
2616  {
2617  e[1] = 1;
2618  }
2619  else if (Dimension() == 3)
2620  {
2621  e[1] = 3;
2622  e[2] = 8;
2623  }
2624
2625  for (int p = 0; p < GetNP(); p++)
2626  {
2627  patchTopo->GetElementEdges(p, edges, orient);
2628
2629  CheckKVDirection(p, kvdir);
2630
2631  for (int d = 0; d < Dimension(); d++)
2632  {
2633  bool flip = false;
2634  if (kvdir[d] == -1) {flip = true;}
2635
2636  // Indices in unique and comprehensive sets of the KnotVector
2637  int iun = edges[e[d]];
2638  int icomp = Dimension()*p+d;
2639
2640  // Check if KnotVectors are of equal order
2641  int o1 = KnotVec(iun)->GetOrder();
2642  int o2 = knotVectorsCompr[icomp]->GetOrder();
2643  int diffo = abs(o1 - o2);
2644
2645  if (diffo)
2646  {
2647  mfem::out << "\norder of knotVectorsCompr " << d << " of patch " << p;
2648  mfem::out << " does not agree with knotVectors " << KnotInd(iun) << "\n";
2649  return false;
2650  }
2651
2652  // Check if Knotvectors have the same knots
2653  if (flip) {knotVectorsCompr[icomp]->Flip();}
2654
2655  KnotVec(iun)->Difference(*(knotVectorsCompr[icomp]), diff);
2656
2657  if (flip) {knotVectorsCompr[icomp]->Flip();}
2658
2659  if (diff.Size() > 0)
2660  {
2661  mfem::out << "\nknotVectorsCompr " << d << " of patch " << p;
2662  mfem::out << " does not agree with knotVectors " << KnotInd(iun) << "\n";
2663  return false;
2664  }
2665  }
2666  }
2667  return true;
2668 }
2669
2671 {
2672  Array<int> edges, orient;
2673
2674  kv.SetSize(Dimension());
2675
2676  if (Dimension() == 1)
2677  {
2678  kv[0] = knotVectorsCompr[Dimension()*p];
2679  }
2680  else if (Dimension() == 2)
2681  {
2682  kv[0] = knotVectorsCompr[Dimension()*p];
2683  kv[1] = knotVectorsCompr[Dimension()*p + 1];
2684  }
2685  else
2686  {
2687  kv[0] = knotVectorsCompr[Dimension()*p];
2688  kv[1] = knotVectorsCompr[Dimension()*p + 1];
2689  kv[2] = knotVectorsCompr[Dimension()*p + 2];
2690  }
2691 }
2692
2694 const
2695 {
2696  Array<int> edges, orient;
2697
2698  kv.SetSize(Dimension());
2699
2700  if (Dimension() == 1)
2701  {
2702  kv[0] = knotVectorsCompr[Dimension()*p];
2703  }
2704  else if (Dimension() == 2)
2705  {
2706  kv[0] = knotVectorsCompr[Dimension()*p];
2707  kv[1] = knotVectorsCompr[Dimension()*p + 1];
2708  }
2709  else
2710  {
2711  kv[0] = knotVectorsCompr[Dimension()*p];
2712  kv[1] = knotVectorsCompr[Dimension()*p + 1];
2713  kv[2] = knotVectorsCompr[Dimension()*p + 2];
2714  }
2715 }
2716
2718 {
2719  Array<int> edges;
2720  Array<int> orient;
2721
2722  kv.SetSize(Dimension() - 1);
2723
2724  if (Dimension() == 2)
2725  {
2726  patchTopo->GetBdrElementEdges(p, edges, orient);
2727  kv[0] = KnotVec(edges[0]);
2728  }
2729  else if (Dimension() == 3)
2730  {
2731  patchTopo->GetBdrElementEdges(p, edges, orient);
2732  kv[0] = KnotVec(edges[0]);
2733  kv[1] = KnotVec(edges[1]);
2734  }
2735 }
2736
2738  int p, Array<const KnotVector *> &kv) const
2739 {
2740  Array<int> edges;
2741  Array<int> orient;
2742
2743  kv.SetSize(Dimension() - 1);
2744
2745  if (Dimension() == 2)
2746  {
2747  patchTopo->GetBdrElementEdges(p, edges, orient);
2748  kv[0] = KnotVec(edges[0]);
2749  }
2750  else if (Dimension() == 3)
2751  {
2752  patchTopo->GetBdrElementEdges(p, edges, orient);
2753  kv[0] = KnotVec(edges[0]);
2754  kv[1] = KnotVec(edges[1]);
2755  }
2756 }
2757
2759 {
2760  MFEM_VERIFY(mOrders.Size() > 0, "");
2761  mOrder = mOrders[0];
2762  for (int i = 1; i < mOrders.Size(); i++)
2763  {
2764  if (mOrders[i] != mOrder)
2765  {
2767  return;
2768  }
2769  }
2770 }
2771
2773 {
2774  mOrders.SetSize(NumOfKnotVectors);
2775  for (int i = 0; i < NumOfKnotVectors; i++)
2776  {
2777  mOrders[i] = knotVectors[i]->GetOrder();
2778  }
2779  SetOrderFromOrders();
2780 }
2781
2783 {
2784  int nv = patchTopo->GetNV();
2785  int ne = patchTopo->GetNEdges();
2786  int nf = patchTopo->GetNFaces();
2787  int np = patchTopo->GetNE();
2788  int meshCounter, spaceCounter, dim = Dimension();
2789
2790  Array<int> edges;
2791  Array<int> orient;
2792
2793  v_meshOffsets.SetSize(nv);
2794  e_meshOffsets.SetSize(ne);
2795  f_meshOffsets.SetSize(nf);
2796  p_meshOffsets.SetSize(np);
2797
2798  v_spaceOffsets.SetSize(nv);
2799  e_spaceOffsets.SetSize(ne);
2800  f_spaceOffsets.SetSize(nf);
2801  p_spaceOffsets.SetSize(np);
2802
2803  // Get vertex offsets
2804  for (meshCounter = 0; meshCounter < nv; meshCounter++)
2805  {
2806  v_meshOffsets[meshCounter] = meshCounter;
2807  v_spaceOffsets[meshCounter] = meshCounter;
2808  }
2809  spaceCounter = meshCounter;
2810
2811  // Get edge offsets
2812  for (int e = 0; e < ne; e++)
2813  {
2814  e_meshOffsets[e] = meshCounter;
2815  e_spaceOffsets[e] = spaceCounter;
2816  meshCounter += KnotVec(e)->GetNE() - 1;
2817  spaceCounter += KnotVec(e)->GetNCP() - 2;
2818  }
2819
2820  // Get face offsets
2821  for (int f = 0; f < nf; f++)
2822  {
2823  f_meshOffsets[f] = meshCounter;
2824  f_spaceOffsets[f] = spaceCounter;
2825
2826  patchTopo->GetFaceEdges(f, edges, orient);
2827
2828  meshCounter +=
2829  (KnotVec(edges[0])->GetNE() - 1) *
2830  (KnotVec(edges[1])->GetNE() - 1);
2831  spaceCounter +=
2832  (KnotVec(edges[0])->GetNCP() - 2) *
2833  (KnotVec(edges[1])->GetNCP() - 2);
2834  }
2835
2836  // Get patch offsets
2837  for (int p = 0; p < np; p++)
2838  {
2839  p_meshOffsets[p] = meshCounter;
2840  p_spaceOffsets[p] = spaceCounter;
2841
2842
2843
2844  if (dim == 1)
2845  {
2846  meshCounter += KnotVec(0)->GetNE() - 1;
2847  spaceCounter += KnotVec(0)->GetNCP() - 2;
2848  }
2849  else if (dim == 2)
2850  {
2851  patchTopo->GetElementEdges(p, edges, orient);
2852  meshCounter +=
2853  (KnotVec(edges[0])->GetNE() - 1) *
2854  (KnotVec(edges[1])->GetNE() - 1);
2855  spaceCounter +=
2856  (KnotVec(edges[0])->GetNCP() - 2) *
2857  (KnotVec(edges[1])->GetNCP() - 2);
2858  }
2859  else
2860  {
2861  patchTopo->GetElementEdges(p, edges, orient);
2862  meshCounter +=
2863  (KnotVec(edges[0])->GetNE() - 1) *
2864  (KnotVec(edges[3])->GetNE() - 1) *
2865  (KnotVec(edges[8])->GetNE() - 1);
2866  spaceCounter +=
2867  (KnotVec(edges[0])->GetNCP() - 2) *
2868  (KnotVec(edges[3])->GetNCP() - 2) *
2869  (KnotVec(edges[8])->GetNCP() - 2);
2870  }
2871  }
2872  NumOfVertices = meshCounter;
2873  NumOfDofs = spaceCounter;
2874 }
2875
2877 {
2878  int dim = Dimension();
2880
2881  NumOfElements = 0;
2882  for (int p = 0; p < GetNP(); p++)
2883  {
2884  GetPatchKnotVectors(p, kv);
2885
2886  int ne = kv[0]->GetNE();
2887  for (int d = 1; d < dim; d++)
2888  {
2889  ne *= kv[d]->GetNE();
2890  }
2891
2892  NumOfElements += ne;
2893  }
2894 }
2895
2897 {
2898  int dim = Dimension() - 1;
2900
2901  NumOfBdrElements = 0;
2902  for (int p = 0; p < GetNBP(); p++)
2903  {
2904  GetBdrPatchKnotVectors(p, kv);
2905
2906  int ne = 1;
2907  for (int d = 0; d < dim; d++)
2908  {
2909  ne *= kv[d]->GetNE();
2910  }
2911
2912  NumOfBdrElements += ne;
2913  }
2914 }
2915
2917 {
2918  elements.SetSize(GetNE());
2919
2920  if (Dimension() == 1)
2921  {
2922  Get1DElementTopo(elements);
2923  }
2924  else if (Dimension() == 2)
2925  {
2926  Get2DElementTopo(elements);
2927  }
2928  else
2929  {
2930  Get3DElementTopo(elements);
2931  }
2932 }
2933
2935 {
2936  int el = 0;
2937  int eg = 0;
2938  int ind[2];
2939  NURBSPatchMap p2g(this);
2940  const KnotVector *kv[1];
2941
2942  for (int p = 0; p < GetNP(); p++)
2943  {
2944  p2g.SetPatchVertexMap(p, kv);
2945  int nx = p2g.nx();
2946
2947  int patch_attr = patchTopo->GetAttribute(p);
2948
2949  for (int i = 0; i < nx; i++)
2950  {
2951  if (activeElem[eg])
2952  {
2953  ind[0] = activeVert[p2g(i)];
2954  ind[1] = activeVert[p2g(i+1)];
2955
2956  elements[el] = new Segment(ind, patch_attr);
2957  el++;
2958  }
2959  eg++;
2960  }
2961  }
2962 }
2963
2965 {
2966  int el = 0;
2967  int eg = 0;
2968  int ind[4];
2969  NURBSPatchMap p2g(this);
2970  const KnotVector *kv[2];
2971
2972  for (int p = 0; p < GetNP(); p++)
2973  {
2974  p2g.SetPatchVertexMap(p, kv);
2975  int nx = p2g.nx();
2976  int ny = p2g.ny();
2977
2978  int patch_attr = patchTopo->GetAttribute(p);
2979
2980  for (int j = 0; j < ny; j++)
2981  {
2982  for (int i = 0; i < nx; i++)
2983  {
2984  if (activeElem[eg])
2985  {
2986  ind[0] = activeVert[p2g(i, j )];
2987  ind[1] = activeVert[p2g(i+1,j )];
2988  ind[2] = activeVert[p2g(i+1,j+1)];
2989  ind[3] = activeVert[p2g(i, j+1)];
2990
2991  elements[el] = new Quadrilateral(ind, patch_attr);
2992  el++;
2993  }
2994  eg++;
2995  }
2996  }
2997  }
2998 }
2999
3001 {
3002  int el = 0;
3003  int eg = 0;
3004  int ind[8];
3005  NURBSPatchMap p2g(this);
3006  const KnotVector *kv[3];
3007
3008  for (int p = 0; p < GetNP(); p++)
3009  {
3010  p2g.SetPatchVertexMap(p, kv);
3011  int nx = p2g.nx();
3012  int ny = p2g.ny();
3013  int nz = p2g.nz();
3014
3015  int patch_attr = patchTopo->GetAttribute(p);
3016
3017  for (int k = 0; k < nz; k++)
3018  {
3019  for (int j = 0; j < ny; j++)
3020  {
3021  for (int i = 0; i < nx; i++)
3022  {
3023  if (activeElem[eg])
3024  {
3025  ind[0] = activeVert[p2g(i, j, k)];
3026  ind[1] = activeVert[p2g(i+1,j, k)];
3027  ind[2] = activeVert[p2g(i+1,j+1,k)];
3028  ind[3] = activeVert[p2g(i, j+1,k)];
3029
3030  ind[4] = activeVert[p2g(i, j, k+1)];
3031  ind[5] = activeVert[p2g(i+1,j, k+1)];
3032  ind[6] = activeVert[p2g(i+1,j+1,k+1)];
3033  ind[7] = activeVert[p2g(i, j+1,k+1)];
3034
3035  elements[el] = new Hexahedron(ind, patch_attr);
3036  el++;
3037  }
3038  eg++;
3039  }
3040  }
3041  }
3042  }
3043 }
3044
3046 {
3047  boundary.SetSize(GetNBE());
3048
3049  if (Dimension() == 1)
3050  {
3051  Get1DBdrElementTopo(boundary);
3052  }
3053  else if (Dimension() == 2)
3054  {
3055  Get2DBdrElementTopo(boundary);
3056  }
3057  else
3058  {
3059  Get3DBdrElementTopo(boundary);
3060  }
3061 }
3062
3064 {
3065  int g_be, l_be;
3066  int ind[2], okv[1];
3067  NURBSPatchMap p2g(this);
3068  const KnotVector *kv[1];
3069
3070  g_be = l_be = 0;
3071  for (int b = 0; b < GetNBP(); b++)
3072  {
3073  p2g.SetBdrPatchVertexMap(b, kv, okv);
3074  int bdr_patch_attr = patchTopo->GetBdrAttribute(b);
3075
3076  if (activeBdrElem[g_be])
3077  {
3078  ind[0] = activeVert[p2g[0]];
3079  boundary[l_be] = new Point(ind, bdr_patch_attr);
3080  l_be++;
3081  }
3082  g_be++;
3083  }
3084 }
3085
3087 {
3088  int g_be, l_be;
3089  int ind[2], okv[1];
3090  NURBSPatchMap p2g(this);
3091  const KnotVector *kv[1];
3092
3093  g_be = l_be = 0;
3094  for (int b = 0; b < GetNBP(); b++)
3095  {
3096  p2g.SetBdrPatchVertexMap(b, kv, okv);
3097  int nx = p2g.nx();
3098
3099  int bdr_patch_attr = patchTopo->GetBdrAttribute(b);
3100
3101  for (int i = 0; i < nx; i++)
3102  {
3103  if (activeBdrElem[g_be])
3104  {
3105  int i_ = (okv[0] >= 0) ? i : (nx - 1 - i);
3106  ind[0] = activeVert[p2g[i_ ]];
3107  ind[1] = activeVert[p2g[i_+1]];
3108
3109  boundary[l_be] = new Segment(ind, bdr_patch_attr);
3110  l_be++;
3111  }
3112  g_be++;
3113  }
3114  }
3115 }
3116
3118 {
3119  int g_be, l_be;
3120  int ind[4], okv[2];
3121  NURBSPatchMap p2g(this);
3122  const KnotVector *kv[2];
3123
3124  g_be = l_be = 0;
3125  for (int b = 0; b < GetNBP(); b++)
3126  {
3127  p2g.SetBdrPatchVertexMap(b, kv, okv);
3128  int nx = p2g.nx();
3129  int ny = p2g.ny();
3130
3131  int bdr_patch_attr = patchTopo->GetBdrAttribute(b);
3132
3133  for (int j = 0; j < ny; j++)
3134  {
3135  int j_ = (okv[1] >= 0) ? j : (ny - 1 - j);
3136  for (int i = 0; i < nx; i++)
3137  {
3138  if (activeBdrElem[g_be])
3139  {
3140  int i_ = (okv[0] >= 0) ? i : (nx - 1 - i);
3141  ind[0] = activeVert[p2g(i_, j_ )];
3142  ind[1] = activeVert[p2g(i_+1,j_ )];
3143  ind[2] = activeVert[p2g(i_+1,j_+1)];
3144  ind[3] = activeVert[p2g(i_, j_+1)];
3145
3146  boundary[l_be] = new Quadrilateral(ind, bdr_patch_attr);
3147  l_be++;
3148  }
3149  g_be++;
3150  }
3151  }
3152  }
3153 }
3154
3156 {
3157  activeDof.SetSize(GetNTotalDof());
3158  activeDof = 0;
3159
3160  if (Dimension() == 1)
3161  {
3162  Generate1DElementDofTable();
3163  }
3164  else if (Dimension() == 2)
3165  {
3166  Generate2DElementDofTable();
3167  }
3168  else
3169  {
3170  Generate3DElementDofTable();
3171  }
3172
3173  SetPatchToElements();
3174
3175  NumOfActiveDofs = 0;
3176  for (int d = 0; d < GetNTotalDof(); d++)
3177  if (activeDof[d])
3178  {
3179  NumOfActiveDofs++;
3180  activeDof[d] = NumOfActiveDofs;
3181  }
3182
3183  int *dof = el_dof->GetJ();
3184  int ndof = el_dof->Size_of_connections();
3185  for (int i = 0; i < ndof; i++)
3186  {
3187  dof[i] = activeDof[dof[i]] - 1;
3188  }
3189 }
3190
3192 {
3193  int el = 0;
3194  int eg = 0;
3195  const KnotVector *kv[2];
3196  NURBSPatchMap p2g(this);
3197
3198  Array<Connection> el_dof_list;
3199  el_to_patch.SetSize(NumOfActiveElems);
3200  el_to_IJK.SetSize(NumOfActiveElems, 2);
3201
3202  for (int p = 0; p < GetNP(); p++)
3203  {
3204  p2g.SetPatchDofMap(p, kv);
3205
3207  const int ord0 = kv[0]->GetOrder();
3208  for (int i = 0; i < kv[0]->GetNKS(); i++)
3209  {
3210  if (kv[0]->isElement(i))
3211  {
3212  if (activeElem[eg])
3213  {
3214  Connection conn(el,0);
3215  for (int ii = 0; ii <= ord0; ii++)
3216  {
3217  conn.to = DofMap(p2g(i+ii));
3218  activeDof[conn.to] = 1;
3219  el_dof_list.Append(conn);
3220  }
3221  el_to_patch[el] = p;
3222  el_to_IJK(el,0) = i;
3223
3224  el++;
3225  }
3226  eg++;
3227  }
3228  }
3229  }
3230  // We must NOT sort el_dof_list in this case.
3231  el_dof = new Table(NumOfActiveElems, el_dof_list);
3232 }
3233
3235 {
3236  int el = 0;
3237  int eg = 0;
3238  const KnotVector *kv[2];
3239  NURBSPatchMap p2g(this);
3240
3241  Array<Connection> el_dof_list;
3242  el_to_patch.SetSize(NumOfActiveElems);
3243  el_to_IJK.SetSize(NumOfActiveElems, 2);
3244
3245  for (int p = 0; p < GetNP(); p++)
3246  {
3247  p2g.SetPatchDofMap(p, kv);
3248
3250  const int ord0 = kv[0]->GetOrder();
3251  const int ord1 = kv[1]->GetOrder();
3252  for (int j = 0; j < kv[1]->GetNKS(); j++)
3253  {
3254  if (kv[1]->isElement(j))
3255  {
3256  for (int i = 0; i < kv[0]->GetNKS(); i++)
3257  {
3258  if (kv[0]->isElement(i))
3259  {
3260  if (activeElem[eg])
3261  {
3262  Connection conn(el,0);
3263  for (int jj = 0; jj <= ord1; jj++)
3264  {
3265  for (int ii = 0; ii <= ord0; ii++)
3266  {
3267  conn.to = DofMap(p2g(i+ii,j+jj));
3268  activeDof[conn.to] = 1;
3269  el_dof_list.Append(conn);
3270  }
3271  }
3272  el_to_patch[el] = p;
3273  el_to_IJK(el,0) = i;
3274  el_to_IJK(el,1) = j;
3275
3276  el++;
3277  }
3278  eg++;
3279  }
3280  }
3281  }
3282  }
3283  }
3284  // We must NOT sort el_dof_list in this case.
3285  el_dof = new Table(NumOfActiveElems, el_dof_list);
3286 }
3287
3289 {
3290  int el = 0;
3291  int eg = 0;
3292  const KnotVector *kv[3];
3293  NURBSPatchMap p2g(this);
3294
3295  Array<Connection> el_dof_list;
3296  el_to_patch.SetSize(NumOfActiveElems);
3297  el_to_IJK.SetSize(NumOfActiveElems, 3);
3298
3299  for (int p = 0; p < GetNP(); p++)
3300  {
3301  p2g.SetPatchDofMap(p, kv);
3302
3304  const int ord0 = kv[0]->GetOrder();
3305  const int ord1 = kv[1]->GetOrder();
3306  const int ord2 = kv[2]->GetOrder();
3307  for (int k = 0; k < kv[2]->GetNKS(); k++)
3308  {
3309  if (kv[2]->isElement(k))
3310  {
3311  for (int j = 0; j < kv[1]->GetNKS(); j++)
3312  {
3313  if (kv[1]->isElement(j))
3314  {
3315  for (int i = 0; i < kv[0]->GetNKS(); i++)
3316  {
3317  if (kv[0]->isElement(i))
3318  {
3319  if (activeElem[eg])
3320  {
3321  Connection conn(el,0);
3322  for (int kk = 0; kk <= ord2; kk++)
3323  {
3324  for (int jj = 0; jj <= ord1; jj++)
3325  {
3326  for (int ii = 0; ii <= ord0; ii++)
3327  {
3328  conn.to = DofMap(p2g(i+ii, j+jj, k+kk));
3329  activeDof[conn.to] = 1;
3330  el_dof_list.Append(conn);
3331  }
3332  }
3333  }
3334
3335  el_to_patch[el] = p;
3336  el_to_IJK(el,0) = i;
3337  el_to_IJK(el,1) = j;
3338  el_to_IJK(el,2) = k;
3339
3340  el++;
3341  }
3342  eg++;
3343  }
3344  }
3345  }
3346  }
3347  }
3348  }
3349  }
3350  // We must NOT sort el_dof_list in this case.
3351  el_dof = new Table(NumOfActiveElems, el_dof_list);
3352 }
3353
3354 void NURBSExtension::GetPatchDofs(const int patch, Array<int> &dofs)
3355 {
3356  const KnotVector *kv[3];
3357  NURBSPatchMap p2g(this);
3358
3359  p2g.SetPatchDofMap(patch, kv);
3360
3361  if (Dimension() == 1)
3362  {
3363  const int nx = kv[0]->GetNCP();
3364  dofs.SetSize(nx);
3365
3366  for (int i=0; i<nx; ++i)
3367  {
3368  dofs[i] = DofMap(p2g(i));
3369  }
3370  }
3371  else if (Dimension() == 2)
3372  {
3373  const int nx = kv[0]->GetNCP();
3374  const int ny = kv[1]->GetNCP();
3375  dofs.SetSize(nx * ny);
3376
3377  for (int j=0; j<ny; ++j)
3378  for (int i=0; i<nx; ++i)
3379  {
3380  dofs[i + (nx * j)] = DofMap(p2g(i, j));
3381  }
3382  }
3383  else if (Dimension() == 3)
3384  {
3385  const int nx = kv[0]->GetNCP();
3386  const int ny = kv[1]->GetNCP();
3387  const int nz = kv[2]->GetNCP();
3388  dofs.SetSize(nx * ny * nz);
3389
3390  for (int k=0; k<nz; ++k)
3391  for (int j=0; j<ny; ++j)
3392  for (int i=0; i<nx; ++i)
3393  {
3394  dofs[i + (nx * (j + (k * ny)))] = DofMap(p2g(i, j, k));
3395  }
3396  }
3397  else
3398  {
3399  MFEM_ABORT("Only 1D/2D/3D supported currently in NURBSExtension::GetPatchDofs");
3400  }
3401 }
3402
3404 {
3405  if (Dimension() == 1)
3406  {
3407  Generate1DBdrElementDofTable();
3408  }
3409  else if (Dimension() == 2)
3410  {
3411  Generate2DBdrElementDofTable();
3412  }
3413  else
3414  {
3415  Generate3DBdrElementDofTable();
3416  }
3417
3418  SetPatchToBdrElements();
3419
3420  int *dof = bel_dof->GetJ();
3421  int ndof = bel_dof->Size_of_connections();
3422  for (int i = 0; i < ndof; i++)
3423  {
3424  dof[i] = activeDof[dof[i]] - 1;
3425  }
3426 }
3427
3429 {
3430  int gbe = 0;
3431  int lbe = 0, okv[1];
3432  const KnotVector *kv[1];
3433  NURBSPatchMap p2g(this);
3434
3435  Array<Connection> bel_dof_list;
3436  bel_to_patch.SetSize(NumOfActiveBdrElems);
3437  bel_to_IJK.SetSize(NumOfActiveBdrElems, 1);
3438
3439  for (int b = 0; b < GetNBP(); b++)
3440  {
3441  p2g.SetBdrPatchDofMap(b, kv, okv);
3443  if (activeBdrElem[gbe])
3444  {
3445  Connection conn(lbe,0);
3446  conn.to = DofMap(p2g[0]);
3447  bel_dof_list.Append(conn);
3448  bel_to_patch[lbe] = b;
3449  bel_to_IJK(lbe,0) = 0;
3450  lbe++;
3451  }
3452  gbe++;
3453  }
3454  // We must NOT sort bel_dof_list in this case.
3455  bel_dof = new Table(NumOfActiveBdrElems, bel_dof_list);
3456 }
3457
3459 {
3460  int gbe = 0;
3461  int lbe = 0, okv[1];
3462  const KnotVector *kv[1];
3463  NURBSPatchMap p2g(this);
3464
3465  Array<Connection> bel_dof_list;
3466  bel_to_patch.SetSize(NumOfActiveBdrElems);
3467  bel_to_IJK.SetSize(NumOfActiveBdrElems, 1);
3468
3469  for (int b = 0; b < GetNBP(); b++)
3470  {
3471  p2g.SetBdrPatchDofMap(b, kv, okv);
3472  const int nx = p2g.nx(); // NCP-1
3474  const int nks0 = kv[0]->GetNKS();
3475  const int ord0 = kv[0]->GetOrder();
3476  for (int i = 0; i < nks0; i++)
3477  {
3478  if (kv[0]->isElement(i))
3479  {
3480  if (activeBdrElem[gbe])
3481  {
3482  Connection conn(lbe,0);
3483  for (int ii = 0; ii <= ord0; ii++)
3484  {
3485  conn.to = DofMap(p2g[(okv[0] >= 0) ? (i+ii) : (nx-i-ii)]);
3486  bel_dof_list.Append(conn);
3487  }
3488  bel_to_patch[lbe] = b;
3489  bel_to_IJK(lbe,0) = (okv[0] >= 0) ? i : (-1-i);
3490  lbe++;
3491  }
3492  gbe++;
3493  }
3494  }
3495  }
3496  // We must NOT sort bel_dof_list in this case.
3497  bel_dof = new Table(NumOfActiveBdrElems, bel_dof_list);
3498 }
3499
3500
3502 {
3503  int gbe = 0;
3504  int lbe = 0, okv[2];
3505  const KnotVector *kv[2];
3506  NURBSPatchMap p2g(this);
3507
3508  Array<Connection> bel_dof_list;
3509  bel_to_patch.SetSize(NumOfActiveBdrElems);
3510  bel_to_IJK.SetSize(NumOfActiveBdrElems, 2);
3511
3512  for (int b = 0; b < GetNBP(); b++)
3513  {
3514  p2g.SetBdrPatchDofMap(b, kv, okv);
3515  const int nx = p2g.nx(); // NCP0-1
3516  const int ny = p2g.ny(); // NCP1-1
3517
3519  const int nks0 = kv[0]->GetNKS();
3520  const int ord0 = kv[0]->GetOrder();
3521  const int nks1 = kv[1]->GetNKS();
3522  const int ord1 = kv[1]->GetOrder();
3523  for (int j = 0; j < nks1; j++)
3524  {
3525  if (kv[1]->isElement(j))
3526  {
3527  for (int i = 0; i < nks0; i++)
3528  {
3529  if (kv[0]->isElement(i))
3530  {
3531  if (activeBdrElem[gbe])
3532  {
3533  Connection conn(lbe,0);
3534  for (int jj = 0; jj <= ord1; jj++)
3535  {
3536  const int jj_ = (okv[1] >= 0) ? (j+jj) : (ny-j-jj);
3537  for (int ii = 0; ii <= ord0; ii++)
3538  {
3539  const int ii_ = (okv[0] >= 0) ? (i+ii) : (nx-i-ii);
3540  conn.to = DofMap(p2g(ii_, jj_));
3541  bel_dof_list.Append(conn);
3542  }
3543  }
3544  bel_to_patch[lbe] = b;
3545  bel_to_IJK(lbe,0) = (okv[0] >= 0) ? i : (-1-i);
3546  bel_to_IJK(lbe,1) = (okv[1] >= 0) ? j : (-1-j);
3547  lbe++;
3548  }
3549  gbe++;
3550  }
3551  }
3552  }
3553  }
3554  }
3555  // We must NOT sort bel_dof_list in this case.
3556  bel_dof = new Table(NumOfActiveBdrElems, bel_dof_list);
3557 }
3558
3560 {
3561  lvert_vert.SetSize(GetNV());
3562  for (int gv = 0; gv < GetGNV(); gv++)
3563  if (activeVert[gv] >= 0)
3564  {
3565  lvert_vert[activeVert[gv]] = gv;
3566  }
3567 }
3568
3570 {
3571  lelem_elem.SetSize(GetNE());
3572  for (int le = 0, ge = 0; ge < GetGNE(); ge++)
3573  if (activeElem[ge])
3574  {
3575  lelem_elem[le++] = ge;
3576  }
3577 }
3578
3579 void NURBSExtension::LoadFE(int i, const FiniteElement *FE) const
3580 {
3581  const NURBSFiniteElement *NURBSFE =
3582  dynamic_cast<const NURBSFiniteElement *>(FE);
3583
3584  if (NURBSFE->GetElement() != i)
3585  {
3586  Array<int> dofs;
3587  NURBSFE->SetIJK(el_to_IJK.GetRow(i));
3588  if (el_to_patch[i] != NURBSFE->GetPatch())
3589  {
3590  GetPatchKnotVectors(el_to_patch[i], NURBSFE->KnotVectors());
3591  NURBSFE->SetPatch(el_to_patch[i]);
3592  NURBSFE->SetOrder();
3593  }
3594  el_dof->GetRow(i, dofs);
3595  weights.GetSubVector(dofs, NURBSFE->Weights());
3596  NURBSFE->SetElement(i);
3597  }
3598 }
3599
3600 void NURBSExtension::LoadBE(int i, const FiniteElement *BE) const
3601 {
3602  if (Dimension() == 1) { return; }
3603
3604  const NURBSFiniteElement *NURBSFE =
3605  dynamic_cast<const NURBSFiniteElement *>(BE);
3606
3607  if (NURBSFE->GetElement() != i)
3608  {
3609  Array<int> dofs;
3610  NURBSFE->SetIJK(bel_to_IJK.GetRow(i));
3611  if (bel_to_patch[i] != NURBSFE->GetPatch())
3612  {
3613  GetBdrPatchKnotVectors(bel_to_patch[i], NURBSFE->KnotVectors());
3614  NURBSFE->SetPatch(bel_to_patch[i]);
3615  NURBSFE->SetOrder();
3616  }
3617  bel_dof->GetRow(i, dofs);
3618  weights.GetSubVector(dofs, NURBSFE->Weights());
3619  NURBSFE->SetElement(i);
3620  }
3621 }
3622
3624 {
3625  delete el_dof;
3626  delete bel_dof;
3627
3628  if (patches.Size() == 0)
3629  {
3630  GetPatchNets(Nodes, Dimension());
3631  }
3632 }
3633
3635 {
3636  if (patches.Size() == 0) { return; }
3637
3638  SetSolutionVector(Nodes, Dimension());
3639  patches.SetSize(0);
3640 }
3641
3643 {
3644  if (patches.Size() == 0)
3645  {
3646  mfem_error("NURBSExtension::SetKnotsFromPatches :"
3647  " No patches available!");
3648  }
3649
3651
3652  for (int p = 0; p < patches.Size(); p++)
3653  {
3654  GetPatchKnotVectors(p, kv);
3655
3656  for (int i = 0; i < kv.Size(); i++)
3657  {
3658  *kv[i] = *patches[p]->GetKV(i);
3659  }
3660  }
3661
3662  UpdateUniqueKV();
3663  SetOrdersFromKnotVectors();
3664
3665  GenerateOffsets();
3666  CountElements();
3667  CountBdrElements();
3668
3669  // all elements must be active
3670  NumOfActiveElems = NumOfElements;
3671  activeElem.SetSize(NumOfElements);
3672  activeElem = true;
3673
3674  GenerateActiveVertices();
3675  InitDofMap();
3676  GenerateElementDofTable();
3677  GenerateActiveBdrElems();
3678  GenerateBdrElementDofTable();
3679
3680  ConnectBoundaries();
3681 }
3682
3683 void NURBSExtension::LoadSolution(std::istream &input, GridFunction &sol) const
3684 {
3685  const FiniteElementSpace *fes = sol.FESpace();
3686  MFEM_VERIFY(fes->GetNURBSext() == this, "");
3687
3688  sol.SetSize(fes->GetVSize());
3689
3690  Array<const KnotVector *> kv(Dimension());
3691  NURBSPatchMap p2g(this);
3692  const int vdim = fes->GetVDim();
3693
3694  for (int p = 0; p < GetNP(); p++)
3695  {
3696  skip_comment_lines(input, '#');
3697
3698  p2g.SetPatchDofMap(p, kv);
3699  const int nx = kv[0]->GetNCP();
3700  const int ny = kv[1]->GetNCP();
3701  const int nz = (kv.Size() == 2) ? 1 : kv[2]->GetNCP();
3702  for (int k = 0; k < nz; k++)
3703  {
3704  for (int j = 0; j < ny; j++)
3705  {
3706  for (int i = 0; i < nx; i++)
3707  {
3708  const int ll = (kv.Size() == 2) ? p2g(i,j) : p2g(i,j,k);
3709  const int l = DofMap(ll);
3710  for (int vd = 0; vd < vdim; vd++)
3711  {
3712  input >> sol(fes->DofToVDof(l,vd));
3713  }
3714  }
3715  }
3716  }
3717  }
3718 }
3719
3720 void NURBSExtension::PrintSolution(const GridFunction &sol, std::ostream &os)
3721 const
3722 {
3723  const FiniteElementSpace *fes = sol.FESpace();
3724  MFEM_VERIFY(fes->GetNURBSext() == this, "");
3725
3726  Array<const KnotVector *> kv(Dimension());
3727  NURBSPatchMap p2g(this);
3728  const int vdim = fes->GetVDim();
3729
3730  for (int p = 0; p < GetNP(); p++)
3731  {
3732  os << "\n# patch " << p << "\n\n";
3733
3734  p2g.SetPatchDofMap(p, kv);
3735  const int nx = kv[0]->GetNCP();
3736  const int ny = kv[1]->GetNCP();
3737  const int nz = (kv.Size() == 2) ? 1 : kv[2]->GetNCP();
3738  for (int k = 0; k < nz; k++)
3739  {
3740  for (int j = 0; j < ny; j++)
3741  {
3742  for (int i = 0; i < nx; i++)
3743  {
3744  const int ll = (kv.Size() == 2) ? p2g(i,j) : p2g(i,j,k);
3745  const int l = DofMap(ll);
3746  os << sol(fes->DofToVDof(l,0));
3747  for (int vd = 1; vd < vdim; vd++)
3748  {
3749  os << ' ' << sol(fes->DofToVDof(l,vd));
3750  }
3751  os << '\n';
3752  }
3753  }
3754  }
3755  }
3756 }
3757
3758 void NURBSExtension::DegreeElevate(int rel_degree, int degree)
3759 {
3760  for (int p = 0; p < patches.Size(); p++)
3761  {
3762  for (int dir = 0; dir < patches[p]->GetNKV(); dir++)
3763  {
3764  int oldd = patches[p]->GetKV(dir)->GetOrder();
3765  int newd = std::min(oldd + rel_degree, degree);
3766  if (newd > oldd)
3767  {
3768  patches[p]->DegreeElevate(dir, newd - oldd);
3769  }
3770  }
3771  }
3772 }
3773
3775 {
3776  for (int p = 0; p < patches.Size(); p++)
3777  {
3778  patches[p]->UniformRefinement();
3779  }
3780 }
3781
3783 {
3784  Array<int> edges;
3785  Array<int> orient;
3786  Array<int> kvdir;
3787
3788  Array<KnotVector *> pkv(Dimension());
3789
3790  for (int p = 0; p < patches.Size(); p++)
3791  {
3792  if (Dimension()==1)
3793  {
3794  pkv[0] = kv[KnotInd(p)];
3795  }
3796  else if (Dimension()==2)
3797  {
3798  patchTopo->GetElementEdges(p, edges, orient);
3799  pkv[0] = kv[KnotInd(edges[0])];
3800  pkv[1] = kv[KnotInd(edges[1])];
3801  }
3802  else if (Dimension()==3)
3803  {
3804  patchTopo->GetElementEdges(p, edges, orient);
3805  pkv[0] = kv[KnotInd(edges[0])];
3806  pkv[1] = kv[KnotInd(edges[3])];
3807  pkv[2] = kv[KnotInd(edges[8])];
3808  }
3809
3810
3811  // Check whether inserted knots should be flipped before inserting.
3812  // Knotvectors are stored in a different array pkvc such that the original
3813  // knots which are inserted are not changed.
3814  // We need those knots for multiple patches so they have to remain original
3815  CheckKVDirection(p, kvdir);
3816
3817  Array<KnotVector *> pkvc(Dimension());
3818  for (int d = 0; d < Dimension(); d++)
3819  {
3820  pkvc[d] = new KnotVector(*(pkv[d]));
3821
3822  if (kvdir[d] == -1)
3823  {
3824  pkvc[d]->Flip();
3825  }
3826  }
3827
3828  patches[p]->KnotInsert(pkvc);
3829  for (int d = 0; d < Dimension(); d++) { delete pkvc[d]; }
3830  }
3831 }
3832
3834 {
3835  Array<int> edges;
3836  Array<int> orient;
3837  Array<int> kvdir;
3838
3839  Array<Vector *> pkv(Dimension());
3840
3841  for (int p = 0; p < patches.Size(); p++)
3842  {
3843  if (Dimension()==1)
3844  {
3845  pkv[0] = kv[KnotInd(p)];
3846  }
3847  else if (Dimension()==2)
3848  {
3849  patchTopo->GetElementEdges(p, edges, orient);
3850  pkv[0] = kv[KnotInd(edges[0])];
3851  pkv[1] = kv[KnotInd(edges[1])];
3852  }
3853  else if (Dimension()==3)
3854  {
3855  patchTopo->GetElementEdges(p, edges, orient);
3856  pkv[0] = kv[KnotInd(edges[0])];
3857  pkv[1] = kv[KnotInd(edges[3])];
3858  pkv[2] = kv[KnotInd(edges[8])];
3859  }
3860
3861
3862  // Check whether inserted knots should be flipped before inserting.
3863  // Knotvectors are stored in a different array pkvc such that the original
3864  // knots which are inserted are not changed.
3865  CheckKVDirection(p, kvdir);
3866
3867  Array<Vector *> pkvc(Dimension());
3868  for (int d = 0; d < Dimension(); d++)
3869  {
3870  pkvc[d] = new Vector(*(pkv[d]));
3871
3872  if (kvdir[d] == -1)
3873  {
3874  // Find flip point, for knotvectors that do not have the domain [0:1]
3875  KnotVector *kva = knotVectorsCompr[Dimension()*p+d];
3876  double apb = (*kva)[0] + (*kva)[kva->Size()-1];
3877
3878  // Flip vector
3879  int size = pkvc[d]->Size();
3880  int ns = ceil(size/2.0);
3881  for (int j = 0; j < ns; j++)
3882  {
3883  double tmp = apb - pkvc[d]->Elem(j);
3884  pkvc[d]->Elem(j) = apb - pkvc[d]->Elem(size-1-j);
3885  pkvc[d]->Elem(size-1-j) = tmp;
3886  }
3887  }
3888  }
3889
3890  patches[p]->KnotInsert(pkvc);
3891
3892  for (int i = 0; i < Dimension(); i++) { delete pkvc[i]; }
3893  }
3894 }
3895
3896 void NURBSExtension::GetPatchNets(const Vector &coords, int vdim)
3897 {
3898  if (Dimension() == 1)
3899  {
3900  Get1DPatchNets(coords, vdim);
3901  }
3902  else if (Dimension() == 2)
3903  {
3904  Get2DPatchNets(coords, vdim);
3905  }
3906  else
3907  {
3908  Get3DPatchNets(coords, vdim);
3909  }
3910 }
3911
3912 void NURBSExtension::Get1DPatchNets(const Vector &coords, int vdim)
3913 {
3915  NURBSPatchMap p2g(this);
3916
3917  patches.SetSize(GetNP());
3918  for (int p = 0; p < GetNP(); p++)
3919  {
3920  p2g.SetPatchDofMap(p, kv);
3921  patches[p] = new NURBSPatch(kv, vdim+1);
3922  NURBSPatch &Patch = *patches[p];
3923
3924  for (int i = 0; i < kv[0]->GetNCP(); i++)
3925  {
3926  const int l = DofMap(p2g(i));
3927  for (int d = 0; d < vdim; d++)
3928  {
3929  Patch(i,d) = coords(l*vdim + d)*weights(l);
3930  }
3931  Patch(i,vdim) = weights(l);
3932  }
3933  }
3934
3935 }
3936
3937 void NURBSExtension::Get2DPatchNets(const Vector &coords, int vdim)
3938 {
3940  NURBSPatchMap p2g(this);
3941
3942  patches.SetSize(GetNP());
3943  for (int p = 0; p < GetNP(); p++)
3944  {
3945  p2g.SetPatchDofMap(p, kv);
3946  patches[p] = new NURBSPatch(kv, vdim+1);
3947  NURBSPatch &Patch = *patches[p];
3948
3949  for (int j = 0; j < kv[1]->GetNCP(); j++)
3950  {
3951  for (int i = 0; i < kv[0]->GetNCP(); i++)
3952  {
3953  const int l = DofMap(p2g(i,j));
3954  for (int d = 0; d < vdim; d++)
3955  {
3956  Patch(i,j,d) = coords(l*vdim + d)*weights(l);
3957  }
3958  Patch(i,j,vdim) = weights(l);
3959  }
3960  }
3961  }
3962 }
3963
3964 void NURBSExtension::Get3DPatchNets(const Vector &coords, int vdim)
3965 {
3967  NURBSPatchMap p2g(this);
3968
3969  patches.SetSize(GetNP());
3970  for (int p = 0; p < GetNP(); p++)
3971  {
3972  p2g.SetPatchDofMap(p, kv);
3973  patches[p] = new NURBSPatch(kv, vdim+1);
3974  NURBSPatch &Patch = *patches[p];
3975
3976  for (int k = 0; k < kv[2]->GetNCP(); k++)
3977  {
3978  for (int j = 0; j < kv[1]->GetNCP(); j++)
3979  {
3980  for (int i = 0; i < kv[0]->GetNCP(); i++)
3981  {
3982  const int l = DofMap(p2g(i,j,k));
3983  for (int d = 0; d < vdim; d++)
3984  {
3985  Patch(i,j,k,d) = coords(l*vdim + d)*weights(l);
3986  }
3987  Patch(i,j,k,vdim) = weights(l);
3988  }
3989  }
3990  }
3991  }
3992 }
3993
3995 {
3996  if (Dimension() == 1)
3997  {
3998  Set1DSolutionVector(coords, vdim);
3999  }
4000  else if (Dimension() == 2)
4001  {
4002  Set2DSolutionVector(coords, vdim);
4003  }
4004  else
4005  {
4006  Set3DSolutionVector(coords, vdim);
4007  }
4008 }
4009
4011 {
4013  NURBSPatchMap p2g(this);
4014
4015  weights.SetSize(GetNDof());
4016  for (int p = 0; p < GetNP(); p++)
4017  {
4018  p2g.SetPatchDofMap(p, kv);
4019  NURBSPatch &patch = *patches[p];
4020  MFEM_ASSERT(vdim+1 == patch.GetNC(), "");
4021
4022  for (int i = 0; i < kv[0]->GetNCP(); i++)
4023  {
4024  const int l = p2g(i);
4025  for (int d = 0; d < vdim; d++)
4026  {
4027  coords(l*vdim + d) = patch(i,d)/patch(i,vdim);
4028  }
4029  weights(l) = patch(i,vdim);
4030  }
4031
4032  delete patches[p];
4033  }
4034 }
4035
4036
4038 {
4040  NURBSPatchMap p2g(this);
4041
4042  weights.SetSize(GetNDof());
4043  for (int p = 0; p < GetNP(); p++)
4044  {
4045  p2g.SetPatchDofMap(p, kv);
4046  NURBSPatch &patch = *patches[p];
4047  MFEM_ASSERT(vdim+1 == patch.GetNC(), "");
4048
4049  for (int j = 0; j < kv[1]->GetNCP(); j++)
4050  {
4051  for (int i = 0; i < kv[0]->GetNCP(); i++)
4052  {
4053  const int l = p2g(i,j);
4054  for (int d = 0; d < vdim; d++)
4055  {
4056  coords(l*vdim + d) = patch(i,j,d)/patch(i,j,vdim);
4057  }
4058  weights(l) = patch(i,j,vdim);
4059  }
4060  }
4061  delete patches[p];
4062  }
4063 }
4064
4066 {
4068  NURBSPatchMap p2g(this);
4069
4070  weights.SetSize(GetNDof());
4071  for (int p = 0; p < GetNP(); p++)
4072  {
4073  p2g.SetPatchDofMap(p, kv);
4074  NURBSPatch &patch = *patches[p];
4075  MFEM_ASSERT(vdim+1 == patch.GetNC(), "");
4076
4077  for (int k = 0; k < kv[2]->GetNCP(); k++)
4078  {
4079  for (int j = 0; j < kv[1]->GetNCP(); j++)
4080  {
4081  for (int i = 0; i < kv[0]->GetNCP(); i++)
4082  {
4083  const int l = p2g(i,j,k);
4084  for (int d = 0; d < vdim; d++)
4085  {
4086  coords(l*vdim + d) = patch(i,j,k,d)/patch(i,j,k,vdim);
4087  }
4088  weights(l) = patch(i,j,k,vdim);
4089  }
4090  }
4091  }
4092  delete patches[p];
4093  }
4094 }
4095
4097 {
4098  MFEM_VERIFY(ijk.Size() == el_to_IJK.NumCols(), "");
4099  el_to_IJK.GetRow(elem, ijk);
4100 }
4101
4103 {
4104  const int np = GetNP();
4105  patch_to_el.resize(np);
4106
4107  for (int e=0; e<el_to_patch.Size(); ++e)
4108  {
4109  patch_to_el[el_to_patch[e]].Append(e);
4110  }
4111 }
4112
4114 {
4115  const int nbp = GetNBP();
4116  patch_to_bel.resize(nbp);
4117
4118  for (int e=0; e<bel_to_patch.Size(); ++e)
4119  {
4120  patch_to_bel[bel_to_patch[e]].Append(e);
4121  }
4122 }
4123
4125 {
4126  MFEM_ASSERT(patch_to_el.size() > 0, "patch_to_el not set");
4127
4128  return patch_to_el[patch];
4129 }
4130
4132 {
4133  MFEM_ASSERT(patch_to_bel.size() > 0, "patch_to_el not set");
4134
4135  return patch_to_bel[patch];
4136 }
4137
4138 #ifdef MFEM_USE_MPI
4140  : NURBSExtension(orig),
4141  partitioning(orig.partitioning ? new int[orig.GetGNE()] : NULL),
4142  gtopo(orig.gtopo),
4143  ldof_group(orig.ldof_group)
4144 {
4145  // Copy the partitioning, if not NULL
4146  if (partitioning)
4147  {
4148  std::memcpy(partitioning, orig.partitioning, orig.GetGNE()*sizeof(int));
4149  }
4150 }
4151
4153  int *part, const Array<bool> &active_bel)
4154  : gtopo(comm)
4155 {
4156  if (parent->NumOfActiveElems < parent->NumOfElements)
4157  {
4158  // SetActive (BuildGroups?) and the way the weights are copied
4159  // do not support this case
4160  mfem_error("ParNURBSExtension::ParNURBSExtension :\n"
4161  " all elements in the parent must be active!");
4162  }
4163
4164  patchTopo = parent->patchTopo;
4165  // steal ownership of patchTopo from the 'parent' NURBS extension
4166  if (!parent->own_topo)
4167  {
4168  mfem_error("ParNURBSExtension::ParNURBSExtension :\n"
4169  " parent does not own the patch topology!");
4170  }
4171  own_topo = 1;
4172  parent->own_topo = 0;
4173
4174  parent->edge_to_knot.Copy(edge_to_knot);
4175
4176  parent->GetOrders().Copy(mOrders);
4177  mOrder = parent->GetOrder();
4178
4179  NumOfKnotVectors = parent->GetNKV();
4180  knotVectors.SetSize(NumOfKnotVectors);
4181  for (int i = 0; i < NumOfKnotVectors; i++)
4182  {
4183  knotVectors[i] = new KnotVector(*parent->GetKnotVector(i));
4184  }
4186
4187  GenerateOffsets();
4188  CountElements();
4189  CountBdrElements();
4190
4191  // copy 'part' to 'partitioning'
4192  partitioning = new int[GetGNE()];
4193  for (int i = 0; i < GetGNE(); i++)
4194  {
4195  partitioning[i] = part[i];
4196  }
4197  SetActive(partitioning, active_bel);
4198
4201  // GenerateActiveBdrElems(); // done by SetActive for now
4203
4204  Table *serial_elem_dof = parent->GetElementDofTable();
4205  BuildGroups(partitioning, *serial_elem_dof);
4206
4207  weights.SetSize(GetNDof());
4208  // copy weights from parent
4209  for (int gel = 0, lel = 0; gel < GetGNE(); gel++)
4210  {
4211  if (activeElem[gel])
4212  {
4213  int ndofs = el_dof->RowSize(lel);
4214  int *ldofs = el_dof->GetRow(lel);
4215  int *gdofs = serial_elem_dof->GetRow(gel);
4216  for (int i = 0; i < ndofs; i++)
4217  {
4218  weights(ldofs[i]) = parent->weights(gdofs[i]);
4219  }
4220  lel++;
4221  }
4222  }
4223 }
4224
4226  const ParNURBSExtension *par_parent)
4227  : gtopo(par_parent->gtopo.GetComm())
4228 {
4229  // steal all data from parent
4230  mOrder = parent->mOrder;
4231  Swap(mOrders, parent->mOrders);
4232
4233  patchTopo = parent->patchTopo;
4234  own_topo = parent->own_topo;
4235  parent->own_topo = 0;
4236
4237  Swap(edge_to_knot, parent->edge_to_knot);
4238
4240  Swap(knotVectors, parent->knotVectors);
4242
4243  NumOfVertices = parent->NumOfVertices;
4244  NumOfElements = parent->NumOfElements;
4246  NumOfDofs = parent->NumOfDofs;
4247
4248  Swap(v_meshOffsets, parent->v_meshOffsets);
4249  Swap(e_meshOffsets, parent->e_meshOffsets);
4250  Swap(f_meshOffsets, parent->f_meshOffsets);
4251  Swap(p_meshOffsets, parent->p_meshOffsets);
4252
4257
4258  Swap(d_to_d, parent->d_to_d);
4259  Swap(master, parent->master);
4260  Swap(slave, parent->slave);
4261
4265  NumOfActiveDofs = parent->NumOfActiveDofs;
4266
4267  Swap(activeVert, parent->activeVert);
4268  Swap(activeElem, parent->activeElem);
4269  Swap(activeBdrElem, parent->activeBdrElem);
4270  Swap(activeDof, parent->activeDof);
4271
4272  el_dof = parent->el_dof;
4273  bel_dof = parent->bel_dof;
4274  parent->el_dof = parent->bel_dof = NULL;
4275
4276  Swap(el_to_patch, parent->el_to_patch);
4277  Swap(bel_to_patch, parent->bel_to_patch);
4278  Swap(el_to_IJK, parent->el_to_IJK);
4279  Swap(bel_to_IJK, parent->bel_to_IJK);
4280
4281  Swap(weights, parent->weights);
4282  MFEM_VERIFY(!parent->HavePatches(), "");
4283
4284  delete parent;
4285
4286  partitioning = NULL;
4287
4288  MFEM_VERIFY(par_parent->partitioning,
4289  "parent ParNURBSExtension has no partitioning!");
4290
4291  // Support for the case when 'parent' is not a local NURBSExtension, i.e.
4292  // NumOfActiveElems is not the same as in 'par_parent'. In that case, we
4293  // assume 'parent' is a global NURBSExtension, i.e. all elements are active.
4294  bool extract_weights = false;
4295  if (NumOfActiveElems != par_parent->NumOfActiveElems)
4296  {
4297  MFEM_ASSERT(NumOfActiveElems == NumOfElements, "internal error");
4298
4299  SetActive(par_parent->partitioning, par_parent->activeBdrElem);
4301  delete el_dof;
4303  el_to_IJK.DeleteAll();
4305  // GenerateActiveBdrElems(); // done by SetActive for now
4306  delete bel_dof;
4310  extract_weights = true;
4311  }
4312
4313  Table *glob_elem_dof = GetGlobalElementDofTable();
4314  BuildGroups(par_parent->partitioning, *glob_elem_dof);
4315  if (extract_weights)
4316  {
4317  Vector glob_weights;
4318  Swap(weights, glob_weights);
4319  weights.SetSize(GetNDof());
4320  // Copy the local 'weights' from the 'glob_weights'.
4321  // Assumption: the local element ids follow the global ordering.
4322  for (int gel = 0, lel = 0; gel < GetGNE(); gel++)
4323  {
4324  if (activeElem[gel])
4325  {
4326  int ndofs = el_dof->RowSize(lel);
4327  int *ldofs = el_dof->GetRow(lel);
4328  int *gdofs = glob_elem_dof->GetRow(gel);
4329  for (int i = 0; i < ndofs; i++)
4330  {
4331  weights(ldofs[i]) = glob_weights(gdofs[i]);
4332  }
4333  lel++;
4334  }
4335  }
4336  }
4337  delete glob_elem_dof;
4338 }
4339
4340 Table *ParNURBSExtension::GetGlobalElementDofTable()
4341 {
4342  if (Dimension() == 1)
4343  {
4344  return Get1DGlobalElementDofTable();
4345  }
4346  else if (Dimension() == 2)
4347  {
4348  return Get2DGlobalElementDofTable();
4349  }
4350  else
4351  {
4352  return Get3DGlobalElementDofTable();
4353  }
4354 }
4355
4356 Table *ParNURBSExtension::Get1DGlobalElementDofTable()
4357 {
4358  int el = 0;
4359  const KnotVector *kv[1];
4360  NURBSPatchMap p2g(this);
4361  Array<Connection> gel_dof_list;
4362
4363  for (int p = 0; p < GetNP(); p++)
4364  {
4365  p2g.SetPatchDofMap(p, kv);
4366
4368  const int ord0 = kv[0]->GetOrder();
4369
4370  for (int i = 0; i < kv[0]->GetNKS(); i++)
4371  {
4372  if (kv[0]->isElement(i))
4373  {
4374  Connection conn(el,0);
4375  for (int ii = 0; ii <= ord0; ii++)
4376  {
4377  conn.to = DofMap(p2g(i+ii));
4378  gel_dof_list.Append(conn);
4379  }
4380  el++;
4381  }
4382  }
4383  }
4384  // We must NOT sort gel_dof_list in this case.
4385  return (new Table(GetGNE(), gel_dof_list));
4386 }
4387
4388 Table *ParNURBSExtension::Get2DGlobalElementDofTable()
4389 {
4390  int el = 0;
4391  const KnotVector *kv[2];
4392  NURBSPatchMap p2g(this);
4393  Array<Connection> gel_dof_list;
4394
4395  for (int p = 0; p < GetNP(); p++)
4396  {
4397  p2g.SetPatchDofMap(p, kv);
4398
4400  const int ord0 = kv[0]->GetOrder();
4401  const int ord1 = kv[1]->GetOrder();
4402  for (int j = 0; j < kv[1]->GetNKS(); j++)
4403  {
4404  if (kv[1]->isElement(j))
4405  {
4406  for (int i = 0; i < kv[0]->GetNKS(); i++)
4407  {
4408  if (kv[0]->isElement(i))
4409  {
4410  Connection conn(el,0);
4411  for (int jj = 0; jj <= ord1; jj++)
4412  {
4413  for (int ii = 0; ii <= ord0; ii++)
4414  {
4415  conn.to = DofMap(p2g(i+ii,j+jj));
4416  gel_dof_list.Append(conn);
4417  }
4418  }
4419  el++;
4420  }
4421  }
4422  }
4423  }
4424  }
4425  // We must NOT sort gel_dof_list in this case.
4426  return (new Table(GetGNE(), gel_dof_list));
4427 }
4428
4429 Table *ParNURBSExtension::Get3DGlobalElementDofTable()
4430 {
4431  int el = 0;
4432  const KnotVector *kv[3];
4433  NURBSPatchMap p2g(this);
4434  Array<Connection> gel_dof_list;
4435
4436  for (int p = 0; p < GetNP(); p++)
4437  {
4438  p2g.SetPatchDofMap(p, kv);
4439
4441  const int ord0 = kv[0]->GetOrder();
4442  const int ord1 = kv[1]->GetOrder();
4443  const int ord2 = kv[2]->GetOrder();
4444  for (int k = 0; k < kv[2]->GetNKS(); k++)
4445  {
4446  if (kv[2]->isElement(k))
4447  {
4448  for (int j = 0; j < kv[1]->GetNKS(); j++)
4449  {
4450  if (kv[1]->isElement(j))
4451  {
4452  for (int i = 0; i < kv[0]->GetNKS(); i++)
4453  {
4454  if (kv[0]->isElement(i))
4455  {
4456  Connection conn(el,0);
4457  for (int kk = 0; kk <= ord2; kk++)
4458  {
4459  for (int jj = 0; jj <= ord1; jj++)
4460  {
4461  for (int ii = 0; ii <= ord0; ii++)
4462  {
4463  conn.to = DofMap(p2g(i+ii,j+jj,k+kk));
4464  gel_dof_list.Append(conn);
4465  }
4466  }
4467  }
4468  el++;
4469  }
4470  }
4471  }
4472  }
4473  }
4474  }
4475  }
4476  // We must NOT sort gel_dof_list in this case.
4477  return (new Table(GetGNE(), gel_dof_list));
4478 }
4479
4480 void ParNURBSExtension::SetActive(const int *partitioning_,
4481  const Array<bool> &active_bel)
4482 {
4484  activeElem = false;
4485  NumOfActiveElems = 0;
4486  const int MyRank = gtopo.MyRank();
4487  for (int i = 0; i < GetGNE(); i++)
4488  if (partitioning_[i] == MyRank)
4489  {
4490  activeElem[i] = true;
4491  NumOfActiveElems++;
4492  }
4493
4494  active_bel.Copy(activeBdrElem);
4495  NumOfActiveBdrElems = 0;
4496  for (int i = 0; i < GetGNBE(); i++)
4497  if (activeBdrElem[i])
4498  {
4500  }
4501 }
4502
4503 void ParNURBSExtension::BuildGroups(const int *partitioning_,
4504  const Table &elem_dof)
4505 {
4506  Table dof_proc;
4507
4508  ListOfIntegerSets groups;
4509  IntegerSet group;
4510
4511  Transpose(elem_dof, dof_proc); // dof_proc is dof_elem
4512  // convert elements to processors
4513  for (int i = 0; i < dof_proc.Size_of_connections(); i++)
4514  {
4515  dof_proc.GetJ()[i] = partitioning_[dof_proc.GetJ()[i]];
4516  }
4517
4518  // the first group is the local one
4519  int MyRank = gtopo.MyRank();
4520  group.Recreate(1, &MyRank);
4521  groups.Insert(group);
4522
4523  int dof = 0;
4525  for (int d = 0; d < GetNTotalDof(); d++)
4526  if (activeDof[d])
4527  {
4528  group.Recreate(dof_proc.RowSize(d), dof_proc.GetRow(d));
4529  ldof_group[dof] = groups.Insert(group);
4530
4531  dof++;
4532  }
4533
4534  gtopo.Create(groups, 1822);
4535 }
4536 #endif // MFEM_USE_MPI
4537
4538
4539 void NURBSPatchMap::GetPatchKnotVectors(int p, const KnotVector *kv[])
4540 {
4541  Ext->patchTopo->GetElementVertices(p, verts);
4542
4543  if (Ext->Dimension() == 1)
4544  {
4545  kv[0] = Ext->knotVectorsCompr[Ext->Dimension()*p];
4546  }
4547  else if (Ext->Dimension() == 2)
4548  {
4549  Ext->patchTopo->GetElementEdges(p, edges, oedge);
4550
4551  kv[0] = Ext->knotVectorsCompr[Ext->Dimension()*p];
4552  kv[1] = Ext->knotVectorsCompr[Ext->Dimension()*p + 1];
4553  }
4554  else if (Ext->Dimension() == 3)
4555  {
4556  Ext->patchTopo->GetElementEdges(p, edges, oedge);
4557  Ext->patchTopo->GetElementFaces(p, faces, oface);
4558
4559  kv[0] = Ext->knotVectorsCompr[Ext->Dimension()*p];
4560  kv[1] = Ext->knotVectorsCompr[Ext->Dimension()*p + 1];
4561  kv[2] = Ext->knotVectorsCompr[Ext->Dimension()*p + 2];
4562  }
4563  opatch = 0;
4564 }
4565
4566 void NURBSPatchMap::GetBdrPatchKnotVectors(int p, const KnotVector *kv[],
4567  int *okv)
4568 {
4569  Ext->patchTopo->GetBdrElementVertices(p, verts);
4570
4571  if (Ext->Dimension() == 2)
4572  {
4573  Ext->patchTopo->GetBdrElementEdges(p, edges, oedge);
4574  kv[0] = Ext->KnotVec(edges[0], oedge[0], &okv[0]);
4575  opatch = oedge[0];
4576  }
4577  else if (Ext->Dimension() == 3)
4578  {
4579  faces.SetSize(1);
4580  Ext->patchTopo->GetBdrElementEdges(p, edges, oedge);
4581  Ext->patchTopo->GetBdrElementFace(p, &faces[0], &opatch);
4582
4583  kv[0] = Ext->KnotVec(edges[0], oedge[0], &okv[0]);
4584  kv[1] = Ext->KnotVec(edges[1], oedge[1], &okv[1]);
4585  }
4586
4587 }
4588
4590 {
4591  GetPatchKnotVectors(p, kv);
4592
4593  I = kv[0]->GetNE() - 1;
4594
4595  for (int i = 0; i < verts.Size(); i++)
4596  {
4597  verts[i] = Ext->v_meshOffsets[verts[i]];
4598  }
4599
4600  if (Ext->Dimension() >= 2)
4601  {
4602  J = kv[1]->GetNE() - 1;
4603  for (int i = 0; i < edges.Size(); i++)
4604  {
4605  edges[i] = Ext->e_meshOffsets[edges[i]];
4606  }
4607  }
4608  if (Ext->Dimension() == 3)
4609  {
4610  K = kv[2]->GetNE() - 1;
4611
4612  for (int i = 0; i < faces.Size(); i++)
4613  {
4614  faces[i] = Ext->f_meshOffsets[faces[i]];
4615  }
4616  }
4617
4618  pOffset = Ext->p_meshOffsets[p];
4619 }
4620
4622 {
4623  GetPatchKnotVectors(p, kv);
4624
4625  I = kv[0]->GetNCP() - 2;
4626
4627  for (int i = 0; i < verts.Size(); i++)
4628  {
4629  verts[i] = Ext->v_spaceOffsets[verts[i]];
4630  }
4631  if (Ext->Dimension() >= 2)
4632  {
4633  J = kv[1]->GetNCP() - 2;
4634  for (int i = 0; i < edges.Size(); i++)
4635  {
4636  edges[i] = Ext->e_spaceOffsets[edges[i]];
4637  }
4638  }
4639  if (Ext->Dimension() == 3)
4640  {
4641  K = kv[2]->GetNCP() - 2;
4642
4643  for (int i = 0; i < faces.Size(); i++)
4644  {
4645  faces[i] = Ext->f_spaceOffsets[faces[i]];
4646  }
4647  }
4648
4649  pOffset = Ext->p_spaceOffsets[p];
4650 }
4651
4653  int *okv)
4654 {
4655  GetBdrPatchKnotVectors(p, kv, okv);
4656
4657  for (int i = 0; i < verts.Size(); i++)
4658  {
4659  verts[i] = Ext->v_meshOffsets[verts[i]];
4660  }
4661
4662  if (Ext->Dimension() == 1)
4663  {
4664  I = 0;
4665  }
4666  else if (Ext->Dimension() == 2)
4667  {
4668  I = kv[0]->GetNE() - 1;
4669  pOffset = Ext->e_meshOffsets[edges[0]];
4670  }
4671  else if (Ext->Dimension() == 3)
4672  {
4673  I = kv[0]->GetNE() - 1;
4674  J = kv[1]->GetNE() - 1;
4675
4676  for (int i = 0; i < edges.Size(); i++)
4677  {
4678  edges[i] = Ext->e_meshOffsets[edges[i]];
4679  }
4680
4681  pOffset = Ext->f_meshOffsets[faces[0]];
4682  }
4683 }
4684
4685 void NURBSPatchMap::SetBdrPatchDofMap(int p, const KnotVector *kv[], int *okv)
4686 {
4687  GetBdrPatchKnotVectors(p, kv, okv);
4688
4689  for (int i = 0; i < verts.Size(); i++)
4690  {
4691  verts[i] = Ext->v_spaceOffsets[verts[i]];
4692  }
4693
4694  if (Ext->Dimension() == 1)
4695  {
4696  I = 0;
4697  }
4698  else if (Ext->Dimension() == 2)
4699  {
4700  I = kv[0]->GetNCP() - 2;
4701  pOffset = Ext->e_spaceOffsets[edges[0]];
4702  }
4703  else if (Ext->Dimension() == 3)
4704  {
4705  I = kv[0]->GetNCP() - 2;
4706  J = kv[1]->GetNCP() - 2;
4707
4708  for (int i = 0; i < edges.Size(); i++)
4709  {
4710  edges[i] = Ext->e_spaceOffsets[edges[i]];
4711  }
4712
4713  pOffset = Ext->f_spaceOffsets[faces[0]];
4714  }
4715 }
4716
4717 }
Array< KnotVector * > knotVectors
Definition: nurbs.hpp:223
Abstract class for all finite elements.
Definition: fe_base.hpp:233
Array< int > f_meshOffsets
Definition: nurbs.hpp:239
void Create(ListOfIntegerSets &groups, int mpitag)
Set up the group topology given the list of sets of shared entities.
Array2D< int > bel_to_IJK
Definition: nurbs.hpp:253
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:605
void SetCoordsFromPatches(Vector &Nodes)
Definition: nurbs.cpp:3634
Array< KnotVector * > knotVectorsCompr
Definition: nurbs.hpp:226
Array< int > activeVert
Definition: nurbs.hpp:214
void ConnectBoundaries2D(int bnd0, int bnd1)
Definition: nurbs.cpp:2080
KnotVector * KnotVec(int edge)
Definition: nurbs.hpp:694
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
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:6298
void DegreeElevate(int dir, int t)
Definition: nurbs.cpp:966
int NumOfElements
Definition: nurbs.hpp:38
NURBSPatch * Interpolate(NURBSPatch &p1, NURBSPatch &p2)
Definition: nurbs.cpp:1398
void Generate2DElementDofTable()
Definition: nurbs.cpp:3234
void MergeGridFunctions(GridFunction *gf_array[], int num_pieces, GridFunction &merged)
Definition: nurbs.cpp:2317
Array< int > ldof_group
Definition: nurbs.hpp:523
bool HavePatches() const
Definition: nurbs.hpp:439
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
Definition: globals.hpp:71
Array< int > master
Definition: nurbs.hpp:233
void PrintSolution(const GridFunction &sol, std::ostream &out) const
Definition: nurbs.cpp:3720
void LoadFE(int i, const FiniteElement *FE) const
Definition: nurbs.cpp:3579
void Generate3DElementDofTable()
Definition: nurbs.cpp:3288
void KnotInsert(int dir, const KnotVector &knot)
Definition: nurbs.cpp:817
void MergeWeights(Mesh *mesh_array[], int num_pieces)
Definition: nurbs.cpp:2292
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:517
void Print(std::ostream &out) const
Definition: nurbs.cpp:126
void CalcDShape(Vector &grad, int i, double xi) const
Definition: nurbs.cpp:188
void UniformRefinement()
Definition: nurbs.cpp:3774
int DofMap(int dof) const
Definition: nurbs.hpp:295
void Get2DBdrElementTopo(Array< Element *> &boundary) const
Definition: nurbs.cpp:3086
void Rotate3D(double normal[], double angle)
Definition: nurbs.cpp:1349
Array< int > e_meshOffsets
Definition: nurbs.hpp:238
int Size() const
Returns the size of the vector.
Definition: vector.hpp:197
void GetElements()
Count the number of elements.
Definition: nurbs.cpp:101
void SetOrdersFromKnotVectors()
Definition: nurbs.cpp:2772
int GetOrder() const
Definition: nurbs.hpp:52
void GetPatchDofs(const int patch, Array< int > &dofs)
Definition: nurbs.cpp:3354
const NURBSExtension * GetNURBSext() const
Definition: fespace.hpp:557
An arbitrary order and dimension NURBS element.
Definition: fe_nurbs.hpp:23
void GetBdrElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of boundary element i.
Definition: mesh.hpp:1297
void Generate1DBdrElementDofTable()
Definition: nurbs.cpp:3428
int Dimension() const
Definition: nurbs.hpp:405
void Generate3DBdrElementDofTable()
Definition: nurbs.cpp:3501
void SetIJK(const int *IJK) const
Definition: fe_nurbs.hpp:50
int GetNE() const
Definition: nurbs.hpp:420
int NumOfControlPoints
Definition: nurbs.hpp:38
Array< int > e_spaceOffsets
Definition: nurbs.hpp:244
void GetBdrElementTopo(Array< Element *> &boundary) const
Definition: nurbs.cpp:3045
void GenerateBdrElementDofTable()
Definition: nurbs.cpp:3403
int MyRank() const
Return the MPI rank within this object&#39;s communicator.
void GetElementIJK(int elem, Array< int > &ijk)
Definition: nurbs.cpp:4096
void LoadSolution(std::istream &input, GridFunction &sol) const
Definition: nurbs.cpp:3683
void GetElementLocalToGlobal(Array< int > &lelem_elem)
Definition: nurbs.cpp:3569
void Get2DElementTopo(Array< Element *> &elements) const
Definition: nurbs.cpp:2964
STL namespace.
void KnotInsert(Array< KnotVector *> &kv)
Definition: nurbs.cpp:3782
void Set1DSolutionVector(Vector &Nodes, int vdim)
Definition: nurbs.cpp:4010
void SetPatchToBdrElements()
Definition: nurbs.cpp:4113
void init(int dim_)
Definition: nurbs.cpp:504
void ConnectBoundaries1D(int bnd0, int bnd1)
Definition: nurbs.cpp:2066
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition: vector.cpp:579
Array< int > edge_to_knot
Definition: nurbs.hpp:221
void GenerateActiveVertices()
Definition: nurbs.cpp:2197
Array< bool > activeElem
Definition: nurbs.hpp:215
void GetVertexLocalToGlobal(Array< int > &lvert_vert)
Definition: nurbs.cpp:3559
void Set2DSolutionVector(Vector &Nodes, int vdim)
Definition: nurbs.cpp:4037
Vector knot
Definition: nurbs.hpp:37
int RowSize(int i) const
Definition: table.hpp:108
void skip_comment_lines(std::istream &is, const char comment_char)
Check if the stream starts with comment_char. If so skip it.
Definition: text.hpp:31
void DeleteAll()
Delete the whole array.
Definition: array.hpp:854
Array< int > el_to_patch
Definition: nurbs.hpp:250
int GetGNBE() const
Definition: nurbs.hpp:421
void SwapDirections(int dir1, int dir2)
Definition: nurbs.cpp:1230
friend class NURBSPatchMap
Definition: nurbs.hpp:202
virtual void SetOrder() const
Update the NURBSFiniteElement according to the currently set knot vectors.
Definition: fe_nurbs.hpp:58
KnotVector()
Create KnotVector.
Definition: nurbs.hpp:42
static const int MaxOrder
Definition: nurbs.hpp:35
Array< int > v_spaceOffsets
Definition: nurbs.hpp:243
std::function< double(const Vector &)> f(double mass_coeff)
Definition: lor_mms.hpp:30
void GetElementTopo(Array< Element *> &elements) const
Definition: nurbs.cpp:2916
void Generate2DBdrElementDofTable()
Definition: nurbs.cpp:3458
Data type hexahedron element.
Definition: hexahedron.hpp:22
void Generate1DElementDofTable()
Definition: nurbs.cpp:3191
void GetBdrElementFace(int i, int *f, int *o) const
Return the index and the orientation of the face of bdr element i. (3D)
Definition: mesh.cpp:6565
void GetBdrPatchKnotVectors(int p, Array< KnotVector *> &kv)
Definition: nurbs.cpp:2717
int GetGNE() const
Definition: nurbs.hpp:419
void ConnectBoundaries3D(int bnd0, int bnd1)
Definition: nurbs.cpp:2126
void GetElementFaces(int i, Array< int > &faces, Array< int > &ori) const
Return the indices and the orientations of all faces of element i.
Definition: mesh.cpp:6514
ParNURBSExtension(const ParNURBSExtension &orig)
Definition: nurbs.cpp:4139
Array< int > d_to_d
Definition: nurbs.hpp:232
void Set3DSolutionVector(Vector &Nodes, int vdim)
Definition: nurbs.cpp:4065
void GenerateActiveBdrElems()
Definition: nurbs.cpp:2270
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:759
int GetOrder() const
If all orders are identical, return that number. Otherwise, return NURBSFECollection::VariableOrder.
Definition: nurbs.hpp:413
void Get3DPatchNets(const Vector &Nodes, int vdim)
Definition: nurbs.cpp:3964
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:154
double b
Definition: lissajous.cpp:42
const Array< int > & GetOrders() const
Definition: nurbs.hpp:410
void Get2DPatchNets(const Vector &Nodes, int vdim)
Definition: nurbs.cpp:3937
int Find(const T &el) const
Return the first index where &#39;el&#39; is found; return -1 if not found.
Definition: array.hpp:818
void GenerateOffsets()
Definition: nurbs.cpp:2782
int GetNCP() const
Definition: nurbs.hpp:51
void UniformRefinement(Vector &newknots) const
Definition: nurbs.cpp:87
Array< int > activeDof
Definition: nurbs.hpp:217
void CalcDnShape(Vector &gradn, int n, int i, double xi) const
Definition: nurbs.cpp:245
void GetElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of element i.
Definition: mesh.hpp:1293
void SetKnotsFromPatches()
Definition: nurbs.cpp:3642
int GetNTotalDof() const
Definition: nurbs.hpp:424
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:6320
Array< const KnotVector * > & KnotVectors() const
Definition: fe_nurbs.hpp:55
Array< int > f_spaceOffsets
Definition: nurbs.hpp:245
Array< bool > activeBdrElem
Definition: nurbs.hpp:216
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:691
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
Definition: table.cpp:187
int GetElement() const
Definition: fe_nurbs.hpp:53
const Array< int > & GetPatchBdrElements(int patch)
Definition: nurbs.cpp:4131
void Get3DElementTopo(Array< Element *> &elements) const
Definition: nurbs.cpp:3000
void Get3DBdrElementTopo(Array< Element *> &boundary) const
Definition: nurbs.cpp:3117
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:702
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for the i&#39;th element. The returned indices are offsets into an ...
Definition: fespace.cpp:281
void SetPatchToElements()
Definition: nurbs.cpp:4102
Array< int > slave
Definition: nurbs.hpp:234
int Size() const
Definition: nurbs.hpp:53
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:555
void Swap(Array< T > &, Array< T > &)
Definition: array.hpp:638
void Difference(const KnotVector &kv, Vector &diff) const
Definition: nurbs.cpp:471
int GetNE() const
Definition: nurbs.hpp:49
void LoadBE(int i, const FiniteElement *BE) const
Definition: nurbs.cpp:3600
void SetElement(int e) const
Definition: fe_nurbs.hpp:54
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:219
int GetNC() const
Definition: nurbs.hpp:152
Array< int > v_meshOffsets
Definition: nurbs.hpp:237
int GetNKS() const
Definition: nurbs.hpp:50
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:66
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:687
void SetBdrPatchVertexMap(int p, const KnotVector *kv[], int *okv)
Definition: nurbs.cpp:4652
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
Definition: table.cpp:413
Array< int > p_meshOffsets
Definition: nurbs.hpp:240
void GetPatchKnotVectors(int p, Array< KnotVector *> &kv)
Definition: nurbs.cpp:2670
KnotVector & operator=(const KnotVector &kv)
Definition: nurbs.cpp:46
void SetPatchVertexMap(int p, const KnotVector *kv[])
Definition: nurbs.cpp:4589
Helper struct for defining a connectivity table, see Table::MakeFromList.
Definition: table.hpp:27
const Array< int > & GetPatchElements(int patch)
Definition: nurbs.cpp:4124
Table * GetElementDofTable()
Definition: nurbs.hpp:442
Array2D< int > el_to_IJK
Definition: nurbs.hpp:252
void SetOrderFromOrders()
Definition: nurbs.cpp:2758
string direction
void Get1DBdrElementTopo(Array< Element *> &boundary) const
Definition: nurbs.cpp:3063
void CountBdrElements()
Definition: nurbs.cpp:2896
void CreateComprehensiveKV()
Definition: nurbs.cpp:2475
void Rotate(double angle, double normal[]=NULL)
Rotate the NURBSPatch.
Definition: nurbs.cpp:1255
void CheckKVDirection(int p, Array< int > &kvdir)
Definition: nurbs.cpp:2407
double a
Definition: lissajous.cpp:41
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:277
void GetPatchNets(const Vector &Nodes, int vdim)
Definition: nurbs.cpp:3896
const KnotVector * GetKnotVector(int i) const
Definition: nurbs.hpp:433
NURBSPatch * Revolve3D(NURBSPatch &patch, double n[], double ang, int times)
Definition: nurbs.cpp:1445
Array< int > p_spaceOffsets
Definition: nurbs.hpp:246
void SetSolutionVector(Vector &Nodes, int vdim)
Definition: nurbs.cpp:3994
void SetBdrPatchDofMap(int p, const KnotVector *kv[], int *okv)
Definition: nurbs.cpp:4685
void ConnectBoundaries()
Definition: nurbs.cpp:2000
int dim
Definition: ex24.cpp:53
void CalcShape(Vector &shape, int i, double xi) const
Definition: nurbs.cpp:161
void GenerateElementDofTable()
Definition: nurbs.cpp:3155
void Copy(Array &copy) const
Create a copy of the internal array to the provided copy.
Definition: array.hpp:864
constexpr int dimension
This example only works in 3D. Kernels for 2D are not implemented.
Definition: hooke.cpp:45
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:709
void Get1DElementTopo(Array< Element *> &elements) const
Definition: nurbs.cpp:2934
int GetNKV() const
Definition: nurbs.hpp:415
Array< int > mOrders
Definition: nurbs.hpp:206
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
RefCoord t[3]
int DofToVDof(int dof, int vd, int ndofs=-1) const
Compute a single vdof corresponding to the index dof and the vector index vd.
Definition: fespace.cpp:251
Vector data type.
Definition: vector.hpp:58
Data type point element.
Definition: point.hpp:22
void DeleteAll()
Delete all dynamically allocated memory, resetting all dimensions to zero.
Definition: array.hpp:444
Vector & Weights() const
Definition: fe_nurbs.hpp:56
bool ConsistentKVSets()
Definition: nurbs.cpp:2603
void Get1DPatchNets(const Vector &Nodes, int vdim)
Definition: nurbs.cpp:3912
void DegreeElevate(int rel_degree, int degree=16)
Definition: nurbs.cpp:3758
void CheckBdrPatches()
Definition: nurbs.cpp:2379
void PrintFunctions(std::ostream &out, int samples=11) const
Definition: nurbs.cpp:133
void SetPatch(int p) const
Definition: fe_nurbs.hpp:52
GroupTopology gtopo
Definition: nurbs.hpp:521
Data type line segment element.
Definition: segment.hpp:22
Array< int > bel_to_patch
Definition: nurbs.hpp:251
void SetPatchDofMap(int p, const KnotVector *kv[])
Definition: nurbs.cpp:4621
int GetNP() const
Definition: nurbs.hpp:406
void ConvertToPatches(const Vector &Nodes)
Definition: nurbs.cpp:3623
int GetNDof() const
Definition: nurbs.hpp:425
int SetLoopDirection(int dir)
Definition: nurbs.cpp:724
KnotVector * DegreeElevate(int t) const
Definition: nurbs.cpp:58