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