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