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