MFEM v4.8.0
Finite element discretization library
Loading...
Searching...
No Matches
fe_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// H1 Finite Element classes utilizing the Bernstein basis
13
14#include "fe_nurbs.hpp"
15#include "../../mesh/nurbs.hpp"
16
17namespace mfem
18{
19
20using namespace std;
21
23{
24 order = kv[0]->GetOrder();
25 dof = order + 1;
26
29}
30
32 Vector &shape) const
33{
34 kv[0]->CalcShape(shape, ijk[0], ip.x);
35
36 real_t sum = 0.0;
37 for (int i = 0; i <= order; i++)
38 {
39 sum += (shape(i) *= weights(i));
40 }
41
42 shape /= sum;
43}
44
46 DenseMatrix &dshape) const
47{
48 Vector grad(dshape.Data(), dof);
49
50 kv[0]->CalcShape (shape_x, ijk[0], ip.x);
51 kv[0]->CalcDShape(grad, ijk[0], ip.x);
52
53 real_t sum = 0.0, dsum = 0.0;
54 for (int i = 0; i <= order; i++)
55 {
56 sum += (shape_x(i) *= weights(i));
57 dsum += ( grad(i) *= weights(i));
58 }
59
60 sum = 1.0/sum;
61 add(sum, grad, -dsum*sum*sum, shape_x, grad);
62}
63
65 DenseMatrix &hessian) const
66{
67 Vector grad(dof);
68 Vector hess(hessian.Data(), dof);
69
70 kv[0]->CalcShape (shape_x, ijk[0], ip.x);
71 kv[0]->CalcDShape(grad, ijk[0], ip.x);
72 kv[0]->CalcD2Shape(hess, ijk[0], ip.x);
73
74 real_t sum = 0.0, dsum = 0.0, d2sum = 0.0;
75 for (int i = 0; i <= order; i++)
76 {
77 sum += (shape_x(i) *= weights(i));
78 dsum += ( grad(i) *= weights(i));
79 d2sum += ( hess(i) *= weights(i));
80 }
81
82 sum = 1.0/sum;
83 add(sum, hess, -2*dsum*sum*sum, grad, hess);
84 add(1.0, hess, (-d2sum + 2*dsum*dsum*sum)*sum*sum, shape_x, hess);
85}
86
87
89{
90 orders[0] = kv[0]->GetOrder();
91 orders[1] = kv[1]->GetOrder();
98
99 order = max(orders[0], orders[1]);
100 dof = (orders[0] + 1)*(orders[1] + 1);
101 u.SetSize(dof);
102 du.SetSize(dof);
104}
105
107 Vector &shape) const
108{
109 kv[0]->CalcShape(shape_x, ijk[0], ip.x);
110 kv[1]->CalcShape(shape_y, ijk[1], ip.y);
111
112 real_t sum = 0.0;
113 for (int o = 0, j = 0; j <= orders[1]; j++)
114 {
115 const real_t sy = shape_y(j);
116 for (int i = 0; i <= orders[0]; i++, o++)
117 {
118 sum += ( shape(o) = shape_x(i)*sy*weights(o) );
119 }
120 }
121
122 shape /= sum;
123}
124
126 DenseMatrix &dshape) const
127{
128 real_t sum, dsum[2];
129
130 kv[0]->CalcShape ( shape_x, ijk[0], ip.x);
131 kv[1]->CalcShape ( shape_y, ijk[1], ip.y);
132
133 kv[0]->CalcDShape(dshape_x, ijk[0], ip.x);
134 kv[1]->CalcDShape(dshape_y, ijk[1], ip.y);
135
136 sum = dsum[0] = dsum[1] = 0.0;
137 for (int o = 0, j = 0; j <= orders[1]; j++)
138 {
139 const real_t sy = shape_y(j), dsy = dshape_y(j);
140 for (int i = 0; i <= orders[0]; i++, o++)
141 {
142 sum += ( u(o) = shape_x(i)*sy*weights(o) );
143
144 dsum[0] += ( dshape(o,0) = dshape_x(i)*sy *weights(o) );
145 dsum[1] += ( dshape(o,1) = shape_x(i)*dsy*weights(o) );
146 }
147 }
148
149 sum = 1.0/sum;
150 dsum[0] *= sum*sum;
151 dsum[1] *= sum*sum;
152
153 for (int o = 0; o < dof; o++)
154 {
155 dshape(o,0) = dshape(o,0)*sum - u(o)*dsum[0];
156 dshape(o,1) = dshape(o,1)*sum - u(o)*dsum[1];
157 }
158}
159
161 DenseMatrix &hessian) const
162{
163 real_t sum, dsum[2], d2sum[3];
164
165 kv[0]->CalcShape ( shape_x, ijk[0], ip.x);
166 kv[1]->CalcShape ( shape_y, ijk[1], ip.y);
167
168 kv[0]->CalcDShape(dshape_x, ijk[0], ip.x);
169 kv[1]->CalcDShape(dshape_y, ijk[1], ip.y);
170
171 kv[0]->CalcD2Shape(d2shape_x, ijk[0], ip.x);
172 kv[1]->CalcD2Shape(d2shape_y, ijk[1], ip.y);
173
174 sum = dsum[0] = dsum[1] = 0.0;
175 d2sum[0] = d2sum[1] = d2sum[2] = 0.0;
176 for (int o = 0, j = 0; j <= orders[1]; j++)
177 {
178 const real_t sy = shape_y(j), dsy = dshape_y(j), d2sy = d2shape_y(j);
179 for (int i = 0; i <= orders[0]; i++, o++)
180 {
181 const real_t sx = shape_x(i), dsx = dshape_x(i), d2sx = d2shape_x(i);
182 sum += ( u(o) = sx*sy*weights(o) );
183
184 dsum[0] += ( du(o,0) = dsx*sy*weights(o) );
185 dsum[1] += ( du(o,1) = sx*dsy*weights(o) );
186
187 d2sum[0] += ( hessian(o,0) = d2sx*sy*weights(o) );
188 d2sum[1] += ( hessian(o,1) = dsx*dsy*weights(o) );
189 d2sum[2] += ( hessian(o,2) = sx*d2sy*weights(o) );
190 }
191 }
192
193 sum = 1.0/sum;
194 dsum[0] *= sum;
195 dsum[1] *= sum;
196
197 d2sum[0] *= sum;
198 d2sum[1] *= sum;
199 d2sum[2] *= sum;
200
201 for (int o = 0; o < dof; o++)
202 {
203 hessian(o,0) = hessian(o,0)*sum
204 - 2*du(o,0)*sum*dsum[0]
205 + u[o]*sum*(2*dsum[0]*dsum[0] - d2sum[0]);
206
207 hessian(o,1) = hessian(o,1)*sum
208 - du(o,0)*sum*dsum[1]
209 - du(o,1)*sum*dsum[0]
210 + u[o]*sum*(2*dsum[0]*dsum[1] - d2sum[1]);
211
212 hessian(o,2) = hessian(o,2)*sum
213 - 2*du(o,1)*sum*dsum[1]
214 + u[o]*sum*(2*dsum[1]*dsum[1] - d2sum[2]);
215 }
216}
217
218
220{
221 orders[0] = kv[0]->GetOrder();
222 orders[1] = kv[1]->GetOrder();
223 orders[2] = kv[2]->GetOrder();
224 shape_x.SetSize(orders[0]+1);
225 shape_y.SetSize(orders[1]+1);
226 shape_z.SetSize(orders[2]+1);
227
228 dshape_x.SetSize(orders[0]+1);
229 dshape_y.SetSize(orders[1]+1);
230 dshape_z.SetSize(orders[2]+1);
231
235
236 order = max(max(orders[0], orders[1]), orders[2]);
237 dof = (orders[0] + 1)*(orders[1] + 1)*(orders[2] + 1);
238 u.SetSize(dof);
239 du.SetSize(dof);
241}
242
244 Vector &shape) const
245{
246 kv[0]->CalcShape(shape_x, ijk[0], ip.x);
247 kv[1]->CalcShape(shape_y, ijk[1], ip.y);
248 kv[2]->CalcShape(shape_z, ijk[2], ip.z);
249
250 real_t sum = 0.0;
251 for (int o = 0, k = 0; k <= orders[2]; k++)
252 {
253 const real_t sz = shape_z(k);
254 for (int j = 0; j <= orders[1]; j++)
255 {
256 const real_t sy_sz = shape_y(j)*sz;
257 for (int i = 0; i <= orders[0]; i++, o++)
258 {
259 sum += ( shape(o) = shape_x(i)*sy_sz*weights(o) );
260 }
261 }
262 }
263
264 shape /= sum;
265}
266
268 DenseMatrix &dshape) const
269{
270 real_t sum, dsum[3];
271
272 kv[0]->CalcShape ( shape_x, ijk[0], ip.x);
273 kv[1]->CalcShape ( shape_y, ijk[1], ip.y);
274 kv[2]->CalcShape ( shape_z, ijk[2], ip.z);
275
276 kv[0]->CalcDShape(dshape_x, ijk[0], ip.x);
277 kv[1]->CalcDShape(dshape_y, ijk[1], ip.y);
278 kv[2]->CalcDShape(dshape_z, ijk[2], ip.z);
279
280 sum = dsum[0] = dsum[1] = dsum[2] = 0.0;
281 for (int o = 0, k = 0; k <= orders[2]; k++)
282 {
283 const real_t sz = shape_z(k), dsz = dshape_z(k);
284 for (int j = 0; j <= orders[1]; j++)
285 {
286 const real_t sy_sz = shape_y(j)* sz;
287 const real_t dsy_sz = dshape_y(j)* sz;
288 const real_t sy_dsz = shape_y(j)*dsz;
289 for (int i = 0; i <= orders[0]; i++, o++)
290 {
291 sum += ( u(o) = shape_x(i)*sy_sz*weights(o) );
292
293 dsum[0] += ( dshape(o,0) = dshape_x(i)* sy_sz *weights(o) );
294 dsum[1] += ( dshape(o,1) = shape_x(i)*dsy_sz *weights(o) );
295 dsum[2] += ( dshape(o,2) = shape_x(i)* sy_dsz*weights(o) );
296 }
297 }
298 }
299
300 sum = 1.0/sum;
301 dsum[0] *= sum*sum;
302 dsum[1] *= sum*sum;
303 dsum[2] *= sum*sum;
304
305 for (int o = 0; o < dof; o++)
306 {
307 dshape(o,0) = dshape(o,0)*sum - u(o)*dsum[0];
308 dshape(o,1) = dshape(o,1)*sum - u(o)*dsum[1];
309 dshape(o,2) = dshape(o,2)*sum - u(o)*dsum[2];
310 }
311}
312
314 DenseMatrix &hessian) const
315{
316 real_t sum, dsum[3], d2sum[6];
317
318 kv[0]->CalcShape ( shape_x, ijk[0], ip.x);
319 kv[1]->CalcShape ( shape_y, ijk[1], ip.y);
320 kv[2]->CalcShape ( shape_z, ijk[2], ip.z);
321
322 kv[0]->CalcDShape(dshape_x, ijk[0], ip.x);
323 kv[1]->CalcDShape(dshape_y, ijk[1], ip.y);
324 kv[2]->CalcDShape(dshape_z, ijk[2], ip.z);
325
326 kv[0]->CalcD2Shape(d2shape_x, ijk[0], ip.x);
327 kv[1]->CalcD2Shape(d2shape_y, ijk[1], ip.y);
328 kv[2]->CalcD2Shape(d2shape_z, ijk[2], ip.z);
329
330 sum = dsum[0] = dsum[1] = dsum[2] = 0.0;
331 d2sum[0] = d2sum[1] = d2sum[2] = d2sum[3] = d2sum[4] = d2sum[5] = 0.0;
332
333 for (int o = 0, k = 0; k <= orders[2]; k++)
334 {
335 const real_t sz = shape_z(k), dsz = dshape_z(k), d2sz = d2shape_z(k);
336 for (int j = 0; j <= orders[1]; j++)
337 {
338 const real_t sy = shape_y(j), dsy = dshape_y(j), d2sy = d2shape_y(j);
339 for (int i = 0; i <= orders[0]; i++, o++)
340 {
341 const real_t sx = shape_x(i), dsx = dshape_x(i), d2sx = d2shape_x(i);
342 sum += ( u(o) = sx*sy*sz*weights(o) );
343
344 dsum[0] += ( du(o,0) = dsx*sy*sz*weights(o) );
345 dsum[1] += ( du(o,1) = sx*dsy*sz*weights(o) );
346 dsum[2] += ( du(o,2) = sx*sy*dsz*weights(o) );
347
348 d2sum[0] += ( hessian(o,0) = d2sx*sy*sz*weights(o) );
349 d2sum[1] += ( hessian(o,1) = dsx*dsy*sz*weights(o) );
350 d2sum[2] += ( hessian(o,2) = dsx*sy*dsz*weights(o) );
351
352 d2sum[3] += ( hessian(o,3) = sx*dsy*dsz*weights(o) );
353
354 d2sum[4] += ( hessian(o,4) = sx*sy*d2sz*weights(o) );
355 d2sum[5] += ( hessian(o,5) = sx*d2sy*sz*weights(o) );
356 }
357 }
358 }
359
360 sum = 1.0/sum;
361 dsum[0] *= sum;
362 dsum[1] *= sum;
363 dsum[2] *= sum;
364
365 d2sum[0] *= sum;
366 d2sum[1] *= sum;
367 d2sum[2] *= sum;
368
369 d2sum[3] *= sum;
370 d2sum[4] *= sum;
371 d2sum[5] *= sum;
372
373 for (int o = 0; o < dof; o++)
374 {
375 hessian(o,0) = hessian(o,0)*sum
376 - 2*du(o,0)*sum*dsum[0]
377 + u[o]*sum*(2*dsum[0]*dsum[0] - d2sum[0]);
378
379 hessian(o,1) = hessian(o,1)*sum
380 - du(o,0)*sum*dsum[1]
381 - du(o,1)*sum*dsum[0]
382 + u[o]*sum*(2*dsum[0]*dsum[1] - d2sum[1]);
383
384 hessian(o,2) = hessian(o,2)*sum
385 - du(o,0)*sum*dsum[2]
386 - du(o,2)*sum*dsum[0]
387 + u[o]*sum*(2*dsum[0]*dsum[2] - d2sum[2]);
388
389 hessian(o,3) = hessian(o,3)*sum
390 - du(o,1)*sum*dsum[2]
391 - du(o,2)*sum*dsum[1]
392 + u[o]*sum*(2*dsum[1]*dsum[2] - d2sum[3]);
393
394 hessian(o,4) = hessian(o,4)*sum
395 - 2*du(o,2)*sum*dsum[2]
396 + u[o]*sum*(2*dsum[2]*dsum[2] - d2sum[4]);
397
398 hessian(o,5) = hessian(o,5)*sum
399 - 2*du(o,1)*sum*dsum[1]
400 + u[o]*sum*(2*dsum[1]*dsum[1] - d2sum[5]);
401 }
402}
403
404
406{
407 orders[0] = kv[0]->GetOrder();
408 orders[1] = kv[1]->GetOrder();
409
410 if (kv1[0]) { delete kv1[0]; }
411 if (kv1[1]) { delete kv1[1]; }
412
413 kv1[0] = kv[0]->DegreeElevate(1);
414 kv1[1] = kv[1]->DegreeElevate(1);
415
416 shape_x.SetSize(orders[0]+1);
417 shape_y.SetSize(orders[1]+1);
418
419 dshape_x.SetSize(orders[0]+1);
420 dshape_y.SetSize(orders[1]+1);
421
424
425 shape1_x.SetSize(orders[0]+2);
426 shape1_y.SetSize(orders[1]+2);
427
430
433
434 order = max(orders[0]+1, orders[1]+1);
435 dof = (orders[0] + 2)*(orders[1] + 1)
436 + (orders[1] + 1)*(orders[1] + 2);
437 u.SetSize(dof);
438 du.SetSize(dof);
440}
441
443 DenseMatrix &shape) const
444{
445 kv[0]->CalcShape(shape_x, ijk[0], ip.x);
446 kv[1]->CalcShape(shape_y, ijk[1], ip.y);
447
448 kv1[0]->CalcShape(shape1_x, ijk[0], ip.x);
449 kv1[1]->CalcShape(shape1_y, ijk[1], ip.y);
450
451 int o = 0;
452 for (int j = 0; j <= orders[1]; j++)
453 {
454 const real_t sy = shape_y(j);
455 for (int i = 0; i <= orders[0]+1; i++, o++)
456 {
457 shape(o,0) = shape1_x(i)*sy;
458 shape(o,1) = 0.0;
459 }
460 }
461
462 for (int j = 0; j <= orders[1]+1; j++)
463 {
464 const real_t sy1 = shape1_y(j);
465 for (int i = 0; i <= orders[0]; i++, o++)
466 {
467 shape(o,0) = 0.0;
468 shape(o,1) = shape_x(i)*sy1;
469 }
470 }
471}
472
474 DenseMatrix &shape) const
475{
476 CalcVShape(Trans.GetIntPoint(), shape);
477 const DenseMatrix & J = Trans.Jacobian();
478 MFEM_ASSERT(J.Width() == 2 && J.Height() == 2,
479 "NURBS_HDiv2DFiniteElement cannot be embedded in "
480 "3 dimensional spaces");
481 for (int i=0; i<dof; i++)
482 {
483 real_t sx = shape(i, 0);
484 real_t sy = shape(i, 1);
485 shape(i, 0) = sx * J(0, 0) + sy * J(0, 1);
486 shape(i, 1) = sx * J(1, 0) + sy * J(1, 1);
487 }
488 shape *= (1.0 / Trans.Weight());
489}
490
492 Vector &divshape) const
493{
494 kv[0]->CalcShape ( shape_x, ijk[0], ip.x);
495 kv[1]->CalcShape ( shape_y, ijk[1], ip.y);
496
497 kv1[0]->CalcDShape(dshape1_x, ijk[0], ip.x);
498 kv1[1]->CalcDShape(dshape1_y, ijk[1], ip.y);
499
500 int o = 0;
501 for (int j = 0; j <= orders[1]; j++)
502 {
503 const real_t sy = shape_y(j);
504 for (int i = 0; i <= orders[0]+1; i++, o++)
505 {
506 divshape(o) = dshape1_x(i)*sy;
507 }
508 }
509
510 for (int j = 0; j <= orders[1]+1; j++)
511 {
512 const real_t dsy1 = dshape1_y(j);
513 for (int i = 0; i <= orders[0]; i++, o++)
514 {
515 divshape(o) = shape_x(i)*dsy1;
516 }
517 }
518}
519
521{
522 if (kv1[0]) { delete kv1[0]; }
523 if (kv1[1]) { delete kv1[1]; }
524}
525
526
528{
529 orders[0] = kv[0]->GetOrder();
530 orders[1] = kv[1]->GetOrder();
531 orders[2] = kv[2]->GetOrder();
532
533 if (kv1[0]) { delete kv1[0]; }
534 if (kv1[1]) { delete kv1[1]; }
535 if (kv1[2]) { delete kv1[2]; }
536
537 kv1[0] = kv[0]->DegreeElevate(1);
538 kv1[1] = kv[1]->DegreeElevate(1);
539 kv1[2] = kv[2]->DegreeElevate(1);
540
541 shape_x.SetSize(orders[0]+1);
542 shape_y.SetSize(orders[1]+1);
543 shape_z.SetSize(orders[2]+1);
544
545 dshape_x.SetSize(orders[0]+1);
546 dshape_y.SetSize(orders[1]+1);
547 dshape_z.SetSize(orders[2]+1);
548
552
553 shape1_x.SetSize(orders[0]+2);
554 shape1_y.SetSize(orders[1]+2);
555 shape1_z.SetSize(orders[2]+2);
556
560
564
565 order = max(orders[0]+1, max( orders[1]+1, orders[2]+1));
566 dof = (orders[0] + 2)*(orders[1] + 1)*(orders[2] + 1) +
567 (orders[0] + 1)*(orders[1] + 2)*(orders[2] + 1) +
568 (orders[0] + 1)*(orders[1] + 1)*(orders[2] + 2);
569 u.SetSize(dof);
570 du.SetSize(dof);
572}
573
575 DenseMatrix &shape) const
576{
577 kv[0]->CalcShape(shape_x, ijk[0], ip.x);
578 kv[1]->CalcShape(shape_y, ijk[1], ip.y);
579 kv[2]->CalcShape(shape_z, ijk[2], ip.z);
580
581 kv1[0]->CalcShape(shape1_x, ijk[0], ip.x);
582 kv1[1]->CalcShape(shape1_y, ijk[1], ip.y);
583 kv1[2]->CalcShape(shape1_z, ijk[2], ip.z);
584
585 shape = 0.0;
586 int o = 0;
587 for (int k = 0; k <= orders[2]; k++)
588 {
589 const real_t sz = shape_z(k);
590 for (int j = 0; j <= orders[1]; j++)
591 {
592 const real_t sy_sz = shape_y(j)*sz;
593 for (int i = 0; i <= orders[0]+1; i++, o++)
594 {
595 shape(o,0) = shape1_x(i)*sy_sz;
596 }
597 }
598 }
599
600 for (int k = 0; k <= orders[2]; k++)
601 {
602 const real_t sz = shape_z(k);
603 for (int j = 0; j <= orders[1]+1; j++)
604 {
605 const real_t sy1_sz = shape1_y(j)*sz;
606 for (int i = 0; i <= orders[0]; i++, o++)
607 {
608 shape(o,1) = shape_x(i)*sy1_sz;
609 }
610 }
611 }
612
613 for (int k = 0; k <= orders[2]+1; k++)
614 {
615 const real_t sz1 = shape1_z(k);
616 for (int j = 0; j <= orders[1]; j++)
617 {
618 const real_t sy_sz1 = shape_y(j)*sz1;
619 for (int i = 0; i <= orders[0]; i++, o++)
620 {
621 shape(o,2) = shape_x(i)*sy_sz1;
622 }
623 }
624 }
625}
626
628 DenseMatrix &shape) const
629{
630 CalcVShape(Trans.GetIntPoint(), shape);
631 const DenseMatrix & J = Trans.Jacobian();
632 MFEM_ASSERT(J.Width() == 3 && J.Height() == 3,
633 "RT_R2D_FiniteElement cannot be embedded in "
634 "3 dimensional spaces");
635 for (int i=0; i<dof; i++)
636 {
637 real_t sx = shape(i, 0);
638 real_t sy = shape(i, 1);
639 real_t sz = shape(i, 2);
640 shape(i, 0) = sx * J(0, 0) + sy * J(0, 1) + sz * J(0, 2);
641 shape(i, 1) = sx * J(1, 0) + sy * J(1, 1) + sz * J(1, 2);
642 shape(i, 2) = sx * J(2, 0) + sy * J(2, 1) + sz * J(2, 2);
643 }
644 shape *= (1.0 / Trans.Weight());
645}
646
648 Vector &divshape) const
649{
650 kv[0]->CalcShape ( shape_x, ijk[0], ip.x);
651 kv[1]->CalcShape ( shape_y, ijk[1], ip.y);
652 kv[2]->CalcShape ( shape_z, ijk[2], ip.z);
653
654 kv1[0]->CalcDShape(dshape1_x, ijk[0], ip.x);
655 kv1[1]->CalcDShape(dshape1_y, ijk[1], ip.y);
656 kv1[2]->CalcDShape(dshape1_z, ijk[2], ip.z);
657
658 int o = 0;
659 for (int k = 0; k <= orders[2]; k++)
660 {
661 const real_t sz = shape_z(k);
662 for (int j = 0; j <= orders[1]; j++)
663 {
664 const real_t sy_sz = shape_y(j)*sz;
665 for (int i = 0; i <= orders[0]+1; i++, o++)
666 {
667 divshape(o) = dshape1_x(i)*sy_sz;
668 }
669 }
670 }
671
672 for (int k = 0; k <= orders[2]; k++)
673 {
674 const real_t sz = shape_z(k);
675 for (int j = 0; j <= orders[1]+1; j++)
676 {
677 const real_t dy1_sz = dshape1_y(j)*sz;
678 for (int i = 0; i <= orders[0]; i++, o++)
679 {
680 divshape(o) = shape_x(i)*dy1_sz;
681 }
682 }
683 }
684
685 for (int k = 0; k <= orders[2]+1; k++)
686 {
687 const real_t dz1 = dshape1_z(k);
688 for (int j = 0; j <= orders[1]; j++)
689 {
690 const real_t sy_dz1 = shape_y(j)*dz1;
691 for (int i = 0; i <= orders[0]; i++, o++)
692 {
693 divshape(o) = shape_x(i)*sy_dz1;
694 }
695 }
696 }
697}
698
700{
701 if (kv1[0]) { delete kv1[0]; }
702 if (kv1[1]) { delete kv1[1]; }
703 if (kv1[2]) { delete kv1[2]; }
704}
705
707{
708 orders[0] = kv[0]->GetOrder();
709 orders[1] = kv[1]->GetOrder();
710
711 if (kv1[0]) { delete kv1[0]; }
712 if (kv1[1]) { delete kv1[1]; }
713
714 kv1[0] = kv[0]->DegreeElevate(1);
715 kv1[1] = kv[1]->DegreeElevate(1);
716
717 shape_x.SetSize(orders[0]+1);
718 shape_y.SetSize(orders[1]+1);
719
720 dshape_x.SetSize(orders[0]+1);
721 dshape_y.SetSize(orders[1]+1);
722
725
726 shape1_x.SetSize(orders[0]+2);
727 shape1_y.SetSize(orders[1]+2);
728
731
734
735 order = max(orders[0]+1, orders[1]+1);
736 dof = (orders[0] + 1)*(orders[1] + 2)
737 + (orders[1] + 2)*(orders[1] + 1);
738 u.SetSize(dof);
739 du.SetSize(dof);
741}
742
744 DenseMatrix &shape) const
745{
746 kv[0]->CalcShape(shape_x, ijk[0], ip.x);
747 kv[1]->CalcShape(shape_y, ijk[1], ip.y);
748
749 kv1[0]->CalcShape(shape1_x, ijk[0], ip.x);
750 kv1[1]->CalcShape(shape1_y, ijk[1], ip.y);
751
752 int o = 0;
753 for (int j = 0; j <= orders[1]+1; j++)
754 {
755 const real_t sy1 = shape1_y(j);
756 for (int i = 0; i <= orders[0]; i++, o++)
757 {
758 shape(o,0) = shape_x(i)*sy1;
759 shape(o,1) = 0.0;
760 }
761 }
762
763 for (int j = 0; j <= orders[1]; j++)
764 {
765 const real_t sy = shape_y(j);
766 for (int i = 0; i <= orders[0]+1; i++, o++)
767 {
768 shape(o,0) = 0.0;
769 shape(o,1) = shape1_x(i)*sy;
770 }
771 }
772}
773
775 DenseMatrix &shape) const
776{
777 CalcVShape(Trans.GetIntPoint(), shape);
778 const DenseMatrix & JI = Trans.InverseJacobian();
779 MFEM_ASSERT(JI.Width() == 2 && JI.Height() == 2,
780 "NURBS_HCurl2DFiniteElement cannot be embedded in "
781 "3 dimensional spaces");
782 for (int i=0; i<dof; i++)
783 {
784 real_t sx = shape(i, 0);
785 real_t sy = shape(i, 1);
786 shape(i, 0) = sx * JI(0, 0) + sy * JI(1, 0);
787 shape(i, 1) = sx * JI(0, 1) + sy * JI(1, 1);
788 }
789}
790
792 DenseMatrix &curl_shape) const
793{
794 kv[0]->CalcShape ( shape_x, ijk[0], ip.x);
795 kv[1]->CalcShape ( shape_y, ijk[1], ip.y);
796
797 kv1[0]->CalcDShape(dshape1_x, ijk[0], ip.x);
798 kv1[1]->CalcDShape(dshape1_y, ijk[1], ip.y);
799
800 int o = 0;
801 for (int j = 0; j <= orders[1]+1; j++)
802 {
803 const real_t dsy1 = dshape1_y(j);
804 for (int i = 0; i <= orders[0]; i++, o++)
805 {
806 curl_shape(o,0) = -shape_x(i)*dsy1;
807 }
808 }
809
810 for (int j = 0; j <= orders[1]; j++)
811 {
812 const real_t sy = shape_y(j);
813 for (int i = 0; i <= orders[0]+1; i++, o++)
814 {
815 curl_shape(o,0) = dshape1_x(i)*sy;
816 }
817 }
818}
819
821{
822 if (kv1[0]) { delete kv1[0]; }
823 if (kv1[1]) { delete kv1[1]; }
824}
825
826
828{
829 orders[0] = kv[0]->GetOrder();
830 orders[1] = kv[1]->GetOrder();
831 orders[2] = kv[2]->GetOrder();
832
833 if (kv1[0]) { delete kv1[0]; }
834 if (kv1[1]) { delete kv1[1]; }
835 if (kv1[2]) { delete kv1[2]; }
836
837 kv1[0] = kv[0]->DegreeElevate(1);
838 kv1[1] = kv[1]->DegreeElevate(1);
839 kv1[2] = kv[2]->DegreeElevate(1);
840
841 shape_x.SetSize(orders[0]+1);
842 shape_y.SetSize(orders[1]+1);
843 shape_z.SetSize(orders[2]+1);
844
845 dshape_x.SetSize(orders[0]+1);
846 dshape_y.SetSize(orders[1]+1);
847 dshape_z.SetSize(orders[2]+1);
848
852
853 shape1_x.SetSize(orders[0]+2);
854 shape1_y.SetSize(orders[1]+2);
855 shape1_z.SetSize(orders[2]+2);
856
860
864
865 order = max(orders[0]+1, max( orders[1]+1, orders[2]+1));
866 dof = (orders[0] + 1)*(orders[1] + 2)*(orders[2] + 2) +
867 (orders[0] + 2)*(orders[1] + 1)*(orders[2] + 2) +
868 (orders[0] + 2)*(orders[1] + 2)*(orders[2] + 1);
869 u.SetSize(dof);
870 du.SetSize(dof);
872}
873
875 DenseMatrix &shape) const
876{
877 kv[0]->CalcShape(shape_x, ijk[0], ip.x);
878 kv[1]->CalcShape(shape_y, ijk[1], ip.y);
879 kv[2]->CalcShape(shape_z, ijk[2], ip.z);
880
881 kv1[0]->CalcShape(shape1_x, ijk[0], ip.x);
882 kv1[1]->CalcShape(shape1_y, ijk[1], ip.y);
883 kv1[2]->CalcShape(shape1_z, ijk[2], ip.z);
884
885 shape = 0.0;
886 int o = 0;
887 for (int k = 0; k <= orders[2]+1; k++)
888 {
889 const real_t sz1 = shape1_z(k);
890 for (int j = 0; j <= orders[1]+1; j++)
891 {
892 const real_t sy1_sz1 = shape1_y(j)*sz1;
893 for (int i = 0; i <= orders[0]; i++, o++)
894 {
895 shape(o,0) = shape_x(i)*sy1_sz1;
896 }
897 }
898 }
899
900 for (int k = 0; k <= orders[2]+1; k++)
901 {
902 const real_t sz1 = shape1_z(k);
903 for (int j = 0; j <= orders[1]; j++)
904 {
905 const real_t sy_sz1 = shape_y(j)*sz1;
906 for (int i = 0; i <= orders[0]+1; i++, o++)
907 {
908 shape(o,1) = shape1_x(i)*sy_sz1;
909 }
910 }
911 }
912
913 for (int k = 0; k <= orders[2]; k++)
914 {
915 const real_t sz = shape_z(k);
916 for (int j = 0; j <= orders[1]+1; j++)
917 {
918 const real_t sy1_sz = shape1_y(j)*sz;
919 for (int i = 0; i <= orders[0]+1; i++, o++)
920 {
921 shape(o,2) = shape1_x(i)*sy1_sz;
922 }
923 }
924 }
925}
926
928 DenseMatrix &shape) const
929{
930 CalcVShape(Trans.GetIntPoint(), shape);
931 const DenseMatrix & JI = Trans.InverseJacobian();
932 MFEM_ASSERT(JI.Width() == 3 && JI.Height() == 3,
933 "NURBS_HCurl3DFiniteElement must be in a"
934 "3 dimensional spaces");
935 for (int i=0; i<dof; i++)
936 {
937 real_t sx = shape(i, 0);
938 real_t sy = shape(i, 1);
939 real_t sz = shape(i, 2);
940 shape(i, 0) = sx * JI(0, 0) + sy * JI(1, 0) + sz * JI(2, 0);
941 shape(i, 1) = sx * JI(0, 1) + sy * JI(1, 1) + sz * JI(2, 1);
942 shape(i, 2) = sx * JI(0, 2) + sy * JI(1, 2) + sz * JI(2, 2);
943 }
944}
945
947 DenseMatrix &curl_shape) const
948{
949 kv[0]->CalcShape ( shape_x, ijk[0], ip.x);
950 kv[1]->CalcShape ( shape_y, ijk[1], ip.y);
951 kv[2]->CalcShape ( shape_z, ijk[2], ip.z);
952
953 kv1[0]->CalcShape(shape1_x, ijk[0], ip.x);
954 kv1[1]->CalcShape(shape1_y, ijk[1], ip.y);
955 kv1[2]->CalcShape(shape1_z, ijk[2], ip.z);
956
957 kv1[0]->CalcDShape(dshape1_x, ijk[0], ip.x);
958 kv1[1]->CalcDShape(dshape1_y, ijk[1], ip.y);
959 kv1[2]->CalcDShape(dshape1_z, ijk[2], ip.z);
960
961 int o = 0;
962 for (int k = 0; k <= orders[2]+1; k++)
963 {
964 const real_t sz1 = shape1_z(k), dsz1 = dshape1_z(k);
965 for (int j = 0; j <= orders[1]+1; j++)
966 {
967 const real_t sy1_dsz1 = shape1_y(j)*dsz1,
968 dsy1_sz1 = dshape1_y(j)*sz1;
969 for (int i = 0; i <= orders[0]; i++, o++)
970 {
971 curl_shape(o,0) = 0.0;
972 curl_shape(o,1) = shape_x(i)*sy1_dsz1;
973 curl_shape(o,2) = -shape_x(i)*dsy1_sz1;
974 }
975 }
976 }
977
978 for (int k = 0; k <= orders[2]+1; k++)
979 {
980 const real_t sz1 = shape1_z(k), dsz1 = dshape1_z(k);
981 for (int j = 0; j <= orders[1]; j++)
982 {
983 const real_t sy_dsz1 = shape_y(j)*dsz1,
984 sy_sz1 = shape_y(j)*sz1;
985 for (int i = 0; i <= orders[0]+1; i++, o++)
986 {
987 curl_shape(o,0) = -shape1_x(i)*sy_dsz1;
988 curl_shape(o,1) = 0.0;
989 curl_shape(o,2) = dshape1_x(i)*sy_sz1;
990 }
991 }
992 }
993
994 for (int k = 0; k <= orders[2]; k++)
995 {
996 const real_t sz = shape_z(k);
997 for (int j = 0; j <= orders[1]+1; j++)
998 {
999 const real_t sy1_sz = shape1_y(j)*sz,
1000 dsy1_sz = dshape1_y(j)*sz;
1001 for (int i = 0; i <= orders[0]+1; i++, o++)
1002 {
1003 curl_shape(o,0) = shape1_x(i)*dsy1_sz;
1004 curl_shape(o,1) = -dshape1_x(i)*sy1_sz;
1005 curl_shape(o,2) = 0.0;
1006 }
1007 }
1008 }
1009}
1010
1012{
1013 if (kv1[0]) { delete kv1[0]; }
1014 if (kv1[1]) { delete kv1[1]; }
1015 if (kv1[2]) { delete kv1[2]; }
1016}
1017
1018}
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
real_t * Data() const
Returns the matrix data array.
Definition densemat.hpp:114
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition densemat.hpp:108
const DenseMatrix & InverseJacobian()
Return the inverse of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition eltrans.hpp:158
const IntegrationPoint & GetIntPoint()
Get a const reference to the currently set integration point. This will return NULL if no integration...
Definition eltrans.hpp:111
real_t Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint....
Definition eltrans.hpp:144
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint,...
Definition eltrans.hpp:132
int dof
Number of degrees of freedom.
Definition fe_base.hpp:253
int orders[Geometry::MaxDim]
Anisotropic orders.
Definition fe_base.hpp:255
int order
Order/degree of the shape functions.
Definition fe_base.hpp:254
Class for integration point with weight.
Definition intrules.hpp:35
void CalcShape(const IntegrationPoint &ip, Vector &shape) const override
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition fe_nurbs.cpp:31
void SetOrder() const override
Definition fe_nurbs.cpp:22
void CalcHessian(const IntegrationPoint &ip, DenseMatrix &hessian) const override
Evaluate the Hessians of all shape functions of a scalar finite element in reference space at the giv...
Definition fe_nurbs.cpp:64
void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const override
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition fe_nurbs.cpp:45
void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const override
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition fe_nurbs.cpp:125
void SetOrder() const override
Definition fe_nurbs.cpp:88
void CalcHessian(const IntegrationPoint &ip, DenseMatrix &hessian) const override
Evaluate the Hessians of all shape functions of a scalar finite element in reference space at the giv...
Definition fe_nurbs.cpp:160
void CalcShape(const IntegrationPoint &ip, Vector &shape) const override
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition fe_nurbs.cpp:106
void SetOrder() const override
Definition fe_nurbs.cpp:219
void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const override
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition fe_nurbs.cpp:267
void CalcHessian(const IntegrationPoint &ip, DenseMatrix &hessian) const override
Evaluate the Hessians of all shape functions of a scalar finite element in reference space at the giv...
Definition fe_nurbs.cpp:313
void CalcShape(const IntegrationPoint &ip, Vector &shape) const override
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition fe_nurbs.cpp:243
Array< const KnotVector * > kv
Definition fe_nurbs.hpp:26
void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const override
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
Definition fe_nurbs.cpp:791
Array< const KnotVector * > kv1
Definition fe_nurbs.hpp:357
void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const override
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Definition fe_nurbs.cpp:743
void SetOrder() const override
Definition fe_nurbs.cpp:706
void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const override
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Definition fe_nurbs.cpp:874
void SetOrder() const override
Definition fe_nurbs.cpp:827
void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const override
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
Definition fe_nurbs.cpp:946
Array< const KnotVector * > kv1
Definition fe_nurbs.hpp:440
void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const override
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Definition fe_nurbs.cpp:442
void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const override
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
Definition fe_nurbs.cpp:491
void SetOrder() const override
Definition fe_nurbs.cpp:405
Array< const KnotVector * > kv1
Definition fe_nurbs.hpp:186
void SetOrder() const override
Definition fe_nurbs.cpp:527
Array< const KnotVector * > kv1
Definition fe_nurbs.hpp:272
void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const override
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
Definition fe_nurbs.cpp:647
void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const override
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Definition fe_nurbs.cpp:574
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition operator.hpp:66
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition operator.hpp:72
Vector data type.
Definition vector.hpp:82
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:558
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition vector.cpp:391
float real_t
Definition config.hpp:43