MFEM v4.10.0
Finite element discretization library
Loading...
Searching...
No Matches
fe_nurbs.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2026, 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
22void NURBS1DFiniteElement::SetOrder() const
23{
24 order = kv[0]->GetOrder();
25 dof = order + 1;
26
27 weights.SetSize(dof);
28 shape_x.SetSize(dof);
29}
30
31void NURBS1DFiniteElement::CalcShape(const IntegrationPoint &ip,
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
45void NURBS1DFiniteElement::CalcDShape(const IntegrationPoint &ip,
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
64void NURBS1DFiniteElement::CalcHessian (const IntegrationPoint &ip,
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
87void NURBS1DFiniteElement::Project(Coefficient &coeff,
88 ElementTransformation &Trans,
89 Vector &dofs) const
90{
91 IntegrationPoint ip;
92
93 for (int i = 0; i <= order; i++)
94 {
95 real_t kx = kv[0]->GetBotella(ijk[0] + i);
96 if (!kv[0]->inSpan(kx, ijk[0]+order)) { continue; }
97 ip.x = kv[0]->GetRefPoint(kx, ijk[0]+order);
98
99 Trans.SetIntPoint(&ip);
100 dofs(i) = coeff.Eval(Trans, ip);
101 }
102}
103
104void NURBS1DFiniteElement::Project(VectorCoefficient &vc,
105 ElementTransformation &Trans,
106 Vector &dofs) const
107{
108 MFEM_ASSERT(dofs.Size() == vc.GetVDim()*dof, "");
109 Vector x(vc.GetVDim());
110 IntegrationPoint ip;
111
112 for (int i = 0; i <= order; i++)
113 {
114 real_t kx = kv[0]->GetBotella(ijk[0] + i);
115 if (!kv[0]->inSpan(kx, ijk[0]+order)) { continue; }
116 ip.x = kv[0]->GetRefPoint(kx, ijk[0]+order);
117
118 Trans.SetIntPoint(&ip);
119 vc.Eval(x, Trans, ip);
120 for (int j = 0; j < x.Size(); j++)
121 {
122 dofs(dof*j+i) = x(j);
123 }
124 }
125}
126
127
128void NURBS2DFiniteElement::SetOrder() const
129{
130 orders[0] = kv[0]->GetOrder();
131 orders[1] = kv[1]->GetOrder();
132 shape_x.SetSize(orders[0]+1);
133 shape_y.SetSize(orders[1]+1);
134 dshape_x.SetSize(orders[0]+1);
135 dshape_y.SetSize(orders[1]+1);
136 d2shape_x.SetSize(orders[0]+1);
137 d2shape_y.SetSize(orders[1]+1);
138
139 order = max(orders[0], orders[1]);
140 dof = (orders[0] + 1)*(orders[1] + 1);
141 u.SetSize(dof);
142 du.SetSize(dof);
143 weights.SetSize(dof);
144}
145
146void NURBS2DFiniteElement::CalcShape(const IntegrationPoint &ip,
147 Vector &shape) const
148{
149 kv[0]->CalcShape(shape_x, ijk[0], ip.x);
150 kv[1]->CalcShape(shape_y, ijk[1], ip.y);
151
152 real_t sum = 0.0;
153 for (int o = 0, j = 0; j <= orders[1]; j++)
154 {
155 const real_t sy = shape_y(j);
156 for (int i = 0; i <= orders[0]; i++, o++)
157 {
158 sum += ( shape(o) = shape_x(i)*sy*weights(o) );
159 }
160 }
161
162 shape /= sum;
163}
164
165void NURBS2DFiniteElement::CalcDShape(const IntegrationPoint &ip,
166 DenseMatrix &dshape) const
167{
168 real_t sum, dsum[2];
169
170 kv[0]->CalcShape ( shape_x, ijk[0], ip.x);
171 kv[1]->CalcShape ( shape_y, ijk[1], ip.y);
172
173 kv[0]->CalcDShape(dshape_x, ijk[0], ip.x);
174 kv[1]->CalcDShape(dshape_y, ijk[1], ip.y);
175
176 sum = dsum[0] = dsum[1] = 0.0;
177 for (int o = 0, j = 0; j <= orders[1]; j++)
178 {
179 const real_t sy = shape_y(j), dsy = dshape_y(j);
180 for (int i = 0; i <= orders[0]; i++, o++)
181 {
182 sum += ( u(o) = shape_x(i)*sy*weights(o) );
183
184 dsum[0] += ( dshape(o,0) = dshape_x(i)*sy *weights(o) );
185 dsum[1] += ( dshape(o,1) = shape_x(i)*dsy*weights(o) );
186 }
187 }
188
189 sum = 1.0/sum;
190 dsum[0] *= sum*sum;
191 dsum[1] *= sum*sum;
192
193 for (int o = 0; o < dof; o++)
194 {
195 dshape(o,0) = dshape(o,0)*sum - u(o)*dsum[0];
196 dshape(o,1) = dshape(o,1)*sum - u(o)*dsum[1];
197 }
198}
199
200void NURBS2DFiniteElement::CalcHessian (const IntegrationPoint &ip,
201 DenseMatrix &hessian) const
202{
203 real_t sum, dsum[2], d2sum[3];
204
205 kv[0]->CalcShape ( shape_x, ijk[0], ip.x);
206 kv[1]->CalcShape ( shape_y, ijk[1], ip.y);
207
208 kv[0]->CalcDShape(dshape_x, ijk[0], ip.x);
209 kv[1]->CalcDShape(dshape_y, ijk[1], ip.y);
210
211 kv[0]->CalcD2Shape(d2shape_x, ijk[0], ip.x);
212 kv[1]->CalcD2Shape(d2shape_y, ijk[1], ip.y);
213
214 sum = dsum[0] = dsum[1] = 0.0;
215 d2sum[0] = d2sum[1] = d2sum[2] = 0.0;
216 for (int o = 0, j = 0; j <= orders[1]; j++)
217 {
218 const real_t sy = shape_y(j), dsy = dshape_y(j), d2sy = d2shape_y(j);
219 for (int i = 0; i <= orders[0]; i++, o++)
220 {
221 const real_t sx = shape_x(i), dsx = dshape_x(i), d2sx = d2shape_x(i);
222 sum += ( u(o) = sx*sy*weights(o) );
223
224 dsum[0] += ( du(o,0) = dsx*sy*weights(o) );
225 dsum[1] += ( du(o,1) = sx*dsy*weights(o) );
226
227 d2sum[0] += ( hessian(o,0) = d2sx*sy*weights(o) );
228 d2sum[1] += ( hessian(o,1) = dsx*dsy*weights(o) );
229 d2sum[2] += ( hessian(o,2) = sx*d2sy*weights(o) );
230 }
231 }
232
233 sum = 1.0/sum;
234 dsum[0] *= sum;
235 dsum[1] *= sum;
236
237 d2sum[0] *= sum;
238 d2sum[1] *= sum;
239 d2sum[2] *= sum;
240
241 for (int o = 0; o < dof; o++)
242 {
243 hessian(o,0) = hessian(o,0)*sum
244 - 2*du(o,0)*sum*dsum[0]
245 + u[o]*sum*(2*dsum[0]*dsum[0] - d2sum[0]);
246
247 hessian(o,1) = hessian(o,1)*sum
248 - du(o,0)*sum*dsum[1]
249 - du(o,1)*sum*dsum[0]
250 + u[o]*sum*(2*dsum[0]*dsum[1] - d2sum[1]);
251
252 hessian(o,2) = hessian(o,2)*sum
253 - 2*du(o,1)*sum*dsum[1]
254 + u[o]*sum*(2*dsum[1]*dsum[1] - d2sum[2]);
255 }
256}
257
258void NURBS2DFiniteElement::Project(Coefficient &coeff,
259 ElementTransformation &Trans,
260 Vector &dofs) const
261{
262 IntegrationPoint ip;
263 for (int o = 0, j = 0; j <= orders[1]; j++)
264 {
265 real_t ky = kv[1]->GetBotella(ijk[1] + j);
266 if (!kv[1]->inSpan(ky, ijk[1]+orders[1]))
267 {
268 o += orders[0] + 1;
269 continue;
270 }
271 ip.y = kv[1]->GetRefPoint(ky, ijk[1]+orders[1]);
272 for (int i = 0; i <= orders[0]; i++, o++)
273 {
274 real_t kx = kv[0]->GetBotella(ijk[0] + i);
275 if (!kv[0]->inSpan(kx, ijk[0]+orders[0])) { continue; }
276 ip.x = kv[0]->GetRefPoint(kx, ijk[0]+orders[0]);
277
278 Trans.SetIntPoint(&ip);
279 dofs(o) = coeff.Eval(Trans, ip);
280 }
281 }
282}
283
284void NURBS2DFiniteElement::Project(VectorCoefficient &vc,
285 ElementTransformation &Trans,
286 Vector &dofs) const
287{
288 MFEM_ASSERT(dofs.Size() == vc.GetVDim()*dof, "");
289 Vector x(vc.GetVDim());
290 IntegrationPoint ip;
291 for (int o = 0, j = 0; j <= orders[1]; j++)
292 {
293 real_t ky = kv[1]->GetBotella(ijk[1] + j);
294 if (!kv[1]->inSpan(ky, ijk[1]+orders[1]))
295 {
296 o += orders[0] + 1;
297 continue;
298 }
299 ip.y = kv[1]->GetRefPoint(ky, ijk[1]+orders[1]);
300 for (int i = 0; i <= orders[0]; i++, o++)
301 {
302 real_t kx = kv[0]->GetBotella(ijk[0] + i);
303 if (!kv[0]->inSpan(kx, ijk[0]+orders[0])) { continue; }
304 ip.x = kv[0]->GetRefPoint(kx, ijk[0]+orders[0]);
305
306 Trans.SetIntPoint(&ip);
307 vc.Eval(x, Trans, ip);
308 for (int v = 0; v < x.Size(); v++)
309 {
310 dofs(dof*v+o) = x(v);
311 }
312 }
313 }
314}
315
316void NURBS3DFiniteElement::SetOrder() const
317{
318 orders[0] = kv[0]->GetOrder();
319 orders[1] = kv[1]->GetOrder();
320 orders[2] = kv[2]->GetOrder();
321 shape_x.SetSize(orders[0]+1);
322 shape_y.SetSize(orders[1]+1);
323 shape_z.SetSize(orders[2]+1);
324
325 dshape_x.SetSize(orders[0]+1);
326 dshape_y.SetSize(orders[1]+1);
327 dshape_z.SetSize(orders[2]+1);
328
329 d2shape_x.SetSize(orders[0]+1);
330 d2shape_y.SetSize(orders[1]+1);
331 d2shape_z.SetSize(orders[2]+1);
332
333 order = max(max(orders[0], orders[1]), orders[2]);
334 dof = (orders[0] + 1)*(orders[1] + 1)*(orders[2] + 1);
335 u.SetSize(dof);
336 du.SetSize(dof);
337 weights.SetSize(dof);
338}
339
340void NURBS3DFiniteElement::CalcShape(const IntegrationPoint &ip,
341 Vector &shape) const
342{
343 kv[0]->CalcShape(shape_x, ijk[0], ip.x);
344 kv[1]->CalcShape(shape_y, ijk[1], ip.y);
345 kv[2]->CalcShape(shape_z, ijk[2], ip.z);
346
347 real_t sum = 0.0;
348 for (int o = 0, k = 0; k <= orders[2]; k++)
349 {
350 const real_t sz = shape_z(k);
351 for (int j = 0; j <= orders[1]; j++)
352 {
353 const real_t sy_sz = shape_y(j)*sz;
354 for (int i = 0; i <= orders[0]; i++, o++)
355 {
356 sum += ( shape(o) = shape_x(i)*sy_sz*weights(o) );
357 }
358 }
359 }
360
361 shape /= sum;
362}
363
364void NURBS3DFiniteElement::CalcDShape(const IntegrationPoint &ip,
365 DenseMatrix &dshape) const
366{
367 real_t sum, dsum[3];
368
369 kv[0]->CalcShape ( shape_x, ijk[0], ip.x);
370 kv[1]->CalcShape ( shape_y, ijk[1], ip.y);
371 kv[2]->CalcShape ( shape_z, ijk[2], ip.z);
372
373 kv[0]->CalcDShape(dshape_x, ijk[0], ip.x);
374 kv[1]->CalcDShape(dshape_y, ijk[1], ip.y);
375 kv[2]->CalcDShape(dshape_z, ijk[2], ip.z);
376
377 sum = dsum[0] = dsum[1] = dsum[2] = 0.0;
378 for (int o = 0, k = 0; k <= orders[2]; k++)
379 {
380 const real_t sz = shape_z(k), dsz = dshape_z(k);
381 for (int j = 0; j <= orders[1]; j++)
382 {
383 const real_t sy_sz = shape_y(j)* sz;
384 const real_t dsy_sz = dshape_y(j)* sz;
385 const real_t sy_dsz = shape_y(j)*dsz;
386 for (int i = 0; i <= orders[0]; i++, o++)
387 {
388 sum += ( u(o) = shape_x(i)*sy_sz*weights(o) );
389
390 dsum[0] += ( dshape(o,0) = dshape_x(i)* sy_sz *weights(o) );
391 dsum[1] += ( dshape(o,1) = shape_x(i)*dsy_sz *weights(o) );
392 dsum[2] += ( dshape(o,2) = shape_x(i)* sy_dsz*weights(o) );
393 }
394 }
395 }
396
397 sum = 1.0/sum;
398 dsum[0] *= sum*sum;
399 dsum[1] *= sum*sum;
400 dsum[2] *= sum*sum;
401
402 for (int o = 0; o < dof; o++)
403 {
404 dshape(o,0) = dshape(o,0)*sum - u(o)*dsum[0];
405 dshape(o,1) = dshape(o,1)*sum - u(o)*dsum[1];
406 dshape(o,2) = dshape(o,2)*sum - u(o)*dsum[2];
407 }
408}
409
410void NURBS3DFiniteElement::CalcHessian (const IntegrationPoint &ip,
411 DenseMatrix &hessian) const
412{
413 real_t sum, dsum[3], d2sum[6];
414
415 kv[0]->CalcShape ( shape_x, ijk[0], ip.x);
416 kv[1]->CalcShape ( shape_y, ijk[1], ip.y);
417 kv[2]->CalcShape ( shape_z, ijk[2], ip.z);
418
419 kv[0]->CalcDShape(dshape_x, ijk[0], ip.x);
420 kv[1]->CalcDShape(dshape_y, ijk[1], ip.y);
421 kv[2]->CalcDShape(dshape_z, ijk[2], ip.z);
422
423 kv[0]->CalcD2Shape(d2shape_x, ijk[0], ip.x);
424 kv[1]->CalcD2Shape(d2shape_y, ijk[1], ip.y);
425 kv[2]->CalcD2Shape(d2shape_z, ijk[2], ip.z);
426
427 sum = dsum[0] = dsum[1] = dsum[2] = 0.0;
428 d2sum[0] = d2sum[1] = d2sum[2] = d2sum[3] = d2sum[4] = d2sum[5] = 0.0;
429
430 for (int o = 0, k = 0; k <= orders[2]; k++)
431 {
432 const real_t sz = shape_z(k), dsz = dshape_z(k), d2sz = d2shape_z(k);
433 for (int j = 0; j <= orders[1]; j++)
434 {
435 const real_t sy = shape_y(j), dsy = dshape_y(j), d2sy = d2shape_y(j);
436 for (int i = 0; i <= orders[0]; i++, o++)
437 {
438 const real_t sx = shape_x(i), dsx = dshape_x(i), d2sx = d2shape_x(i);
439 sum += ( u(o) = sx*sy*sz*weights(o) );
440
441 dsum[0] += ( du(o,0) = dsx*sy*sz*weights(o) );
442 dsum[1] += ( du(o,1) = sx*dsy*sz*weights(o) );
443 dsum[2] += ( du(o,2) = sx*sy*dsz*weights(o) );
444
445 d2sum[0] += ( hessian(o,0) = d2sx*sy*sz*weights(o) );
446 d2sum[1] += ( hessian(o,1) = dsx*dsy*sz*weights(o) );
447 d2sum[2] += ( hessian(o,2) = dsx*sy*dsz*weights(o) );
448 d2sum[3] += ( hessian(o,3) = sx*d2sy*sz*weights(o) );
449 d2sum[4] += ( hessian(o,4) = sx*dsy*dsz*weights(o) );
450 d2sum[5] += ( hessian(o,5) = sx*sy*d2sz*weights(o) );
451
452 }
453 }
454 }
455
456 sum = 1.0/sum;
457 dsum[0] *= sum;
458 dsum[1] *= sum;
459 dsum[2] *= sum;
460
461 d2sum[0] *= sum;
462 d2sum[1] *= sum;
463 d2sum[2] *= sum;
464
465 d2sum[3] *= sum;
466 d2sum[4] *= sum;
467 d2sum[5] *= sum;
468
469 for (int o = 0; o < dof; o++)
470 {
471 hessian(o,0) = hessian(o,0)*sum
472 - 2*du(o,0)*sum*dsum[0]
473 + u[o]*sum*(2*dsum[0]*dsum[0] - d2sum[0]);
474
475 hessian(o,1) = hessian(o,1)*sum
476 - du(o,0)*sum*dsum[1]
477 - du(o,1)*sum*dsum[0]
478 + u[o]*sum*(2*dsum[0]*dsum[1] - d2sum[1]);
479
480 hessian(o,2) = hessian(o,2)*sum
481 - du(o,0)*sum*dsum[2]
482 - du(o,2)*sum*dsum[0]
483 + u[o]*sum*(2*dsum[0]*dsum[2] - d2sum[2]);
484
485 hessian(o,3) = hessian(o,3)*sum
486 - du(o,1)*sum*dsum[2]
487 - du(o,2)*sum*dsum[1]
488 + u[o]*sum*(2*dsum[1]*dsum[2] - d2sum[3]);
489
490 hessian(o,4) = hessian(o,4)*sum
491 - 2*du(o,2)*sum*dsum[2]
492 + u[o]*sum*(2*dsum[2]*dsum[2] - d2sum[4]);
493
494 hessian(o,5) = hessian(o,5)*sum
495 - 2*du(o,1)*sum*dsum[1]
496 + u[o]*sum*(2*dsum[1]*dsum[1] - d2sum[5]);
497 }
498}
499
500void NURBS3DFiniteElement::Project(Coefficient &coeff,
501 ElementTransformation &Trans,
502 Vector &dofs) const
503{
504 IntegrationPoint ip;
505
506 for (int o = 0, k = 0; k <= orders[2]; k++)
507 {
508 real_t kz = kv[2]->GetBotella(ijk[2] + k);
509 if (!kv[2]->inSpan(kz, ijk[2]+orders[2]))
510 {
511 o += (orders[0] + 1)*(orders[1] + 1);
512 continue;
513 }
514 ip.z = kv[2]->GetRefPoint(kz, ijk[2]+orders[2]);
515 for (int j = 0; j <= orders[1]; j++)
516 {
517 real_t ky = kv[1]->GetBotella(ijk[1] + j);
518 if (!kv[1]->inSpan(ky, ijk[1]+orders[1]))
519 {
520 o += orders[0] + 1;
521 continue;
522 }
523 ip.y = kv[1]->GetRefPoint(ky, ijk[1]+orders[1]);
524 for (int i = 0; i <= orders[0]; i++, o++)
525 {
526 real_t kx = kv[0]->GetBotella(ijk[0] + i);
527 if (!kv[0]->inSpan(kx, ijk[0]+orders[0])) { continue; }
528 ip.x = kv[0]->GetRefPoint(kx, ijk[0]+orders[0]);
529
530 Trans.SetIntPoint(&ip);
531 dofs(o) = coeff.Eval(Trans, ip);
532 }
533 }
534 }
535}
536
537void NURBS3DFiniteElement::Project(VectorCoefficient &vc,
538 ElementTransformation &Trans,
539 Vector &dofs) const
540{
541 MFEM_ASSERT(dofs.Size() == vc.GetVDim()*dof, "");
542 Vector x(vc.GetVDim());
543 IntegrationPoint ip;
544
545 for (int o = 0, k = 0; k <= orders[2]; k++)
546 {
547 real_t kz = kv[2]->GetBotella(ijk[2] + k);
548 if (!kv[2]->inSpan(kz, ijk[2]+orders[2]))
549 {
550 o += (orders[0] + 1)*(orders[1] + 1);
551 continue;
552 }
553 ip.z = kv[2]->GetRefPoint(kz, ijk[2]+orders[2]);
554 for (int j = 0; j <= orders[1]; j++)
555 {
556 real_t ky = kv[1]->GetBotella(ijk[1] + j);
557 if (!kv[1]->inSpan(ky, ijk[1]+orders[1]))
558 {
559 o += orders[0] + 1;
560 continue;
561 }
562 ip.y = kv[1]->GetRefPoint(ky, ijk[1]+orders[1]);
563 for (int i = 0; i <= orders[0]; i++, o++)
564 {
565 real_t kx = kv[0]->GetBotella(ijk[0] + i);
566 if (!kv[0]->inSpan(kx, ijk[0]+orders[0])) { continue; }
567 ip.x = kv[0]->GetRefPoint(kx, ijk[0]+orders[0]);
568
569 Trans.SetIntPoint(&ip);
570 vc.Eval(x, Trans, ip);
571 for (int v = 0; v < x.Size(); v++)
572 {
573 dofs(dof*v+o) = x(v);
574 }
575 }
576 }
577 }
578}
579
580void NURBS_HDiv2DFiniteElement::SetOrder() const
581{
582 orders[0] = kv[0]->GetOrder();
583 orders[1] = kv[1]->GetOrder();
584
585 if (kv1[0]) { delete kv1[0]; }
586 if (kv1[1]) { delete kv1[1]; }
587
588 kv1[0] = kv[0]->DegreeElevate(1);
589 kv1[1] = kv[1]->DegreeElevate(1);
590
591 shape_x.SetSize(orders[0]+1);
592 shape_y.SetSize(orders[1]+1);
593
594 dshape_x.SetSize(orders[0]+1);
595 dshape_y.SetSize(orders[1]+1);
596
597 d2shape_x.SetSize(orders[0]+1);
598 d2shape_y.SetSize(orders[1]+1);
599
600 shape1_x.SetSize(orders[0]+2);
601 shape1_y.SetSize(orders[1]+2);
602
603 dshape1_x.SetSize(orders[0]+2);
604 dshape1_y.SetSize(orders[1]+2);
605
606 d2shape1_x.SetSize(orders[0]+2);
607 d2shape1_y.SetSize(orders[1]+2);
608
609 order = max(orders[0]+1, orders[1]+1);
610 dof = (orders[0] + 2)*(orders[1] + 1)
611 + (orders[1] + 1)*(orders[1] + 2);
612 u.SetSize(dof);
613 du.SetSize(dof);
614 weights.SetSize(dof);
615}
616
617void NURBS_HDiv2DFiniteElement::CalcVShape(const IntegrationPoint &ip,
618 DenseMatrix &shape) const
619{
620 kv[0]->CalcShape(shape_x, ijk[0], ip.x);
621 kv[1]->CalcShape(shape_y, ijk[1], ip.y);
622
623 kv1[0]->CalcShape(shape1_x, ijk[0], ip.x);
624 kv1[1]->CalcShape(shape1_y, ijk[1], ip.y);
625
626 int o = 0;
627 for (int j = 0; j <= orders[1]; j++)
628 {
629 const real_t sy = shape_y(j);
630 for (int i = 0; i <= orders[0]+1; i++, o++)
631 {
632 shape(o,0) = shape1_x(i)*sy;
633 shape(o,1) = 0.0;
634 }
635 }
636
637 for (int j = 0; j <= orders[1]+1; j++)
638 {
639 const real_t sy1 = shape1_y(j);
640 for (int i = 0; i <= orders[0]; i++, o++)
641 {
642 shape(o,0) = 0.0;
643 shape(o,1) = shape_x(i)*sy1;
644 }
645 }
646}
647
648void NURBS_HDiv2DFiniteElement::CalcVShape(ElementTransformation &Trans,
649 DenseMatrix &shape) const
650{
651 CalcVShape(Trans.GetIntPoint(), shape);
652 const DenseMatrix & J = Trans.Jacobian();
653 MFEM_ASSERT(J.Width() == 2 && J.Height() == 2,
654 "NURBS_HDiv2DFiniteElement cannot be embedded in "
655 "3 dimensional spaces");
656 for (int i=0; i<dof; i++)
657 {
658 real_t sx = shape(i, 0);
659 real_t sy = shape(i, 1);
660 shape(i, 0) = sx * J(0, 0) + sy * J(0, 1);
661 shape(i, 1) = sx * J(1, 0) + sy * J(1, 1);
662 }
663 shape *= (1.0 / Trans.Weight());
664}
665
666void NURBS_HDiv2DFiniteElement::CalcDivShape(const IntegrationPoint &ip,
667 Vector &divshape) const
668{
669 kv[0]->CalcShape ( shape_x, ijk[0], ip.x);
670 kv[1]->CalcShape ( shape_y, ijk[1], ip.y);
671
672 kv1[0]->CalcDShape(dshape1_x, ijk[0], ip.x);
673 kv1[1]->CalcDShape(dshape1_y, ijk[1], ip.y);
674
675 int o = 0;
676 for (int j = 0; j <= orders[1]; j++)
677 {
678 const real_t sy = shape_y(j);
679 for (int i = 0; i <= orders[0]+1; i++, o++)
680 {
681 divshape(o) = dshape1_x(i)*sy;
682 }
683 }
684
685 for (int j = 0; j <= orders[1]+1; j++)
686 {
687 const real_t dsy1 = dshape1_y(j);
688 for (int i = 0; i <= orders[0]; i++, o++)
689 {
690 divshape(o) = shape_x(i)*dsy1;
691 }
692 }
693}
694
695void NURBS_HDiv2DFiniteElement::Project(VectorCoefficient &vc,
696 ElementTransformation &Trans,
697 Vector &dofs) const
698{
699 MFEM_ASSERT(dofs.Size() == dof, "");
700 MFEM_ASSERT(vc.GetVDim() == 2, "");
701 Vector x(2), mx(2);
702 IntegrationPoint ip;
703 int o = 0;
704
705 for (int j = 0; j <= orders[1]; j++)
706 {
707 real_t ky = kv[1]->GetBotella(ijk[1] + j);
708 if (!kv[1]->inSpan(ky, ijk[1]+orders[1]))
709 {
710 o += orders[0] + 2;
711 continue;
712 }
713 ip.y = kv[1]->GetRefPoint(ky, ijk[1]+orders[1]);
714 for (int i = 0; i <= orders[0]+1; i++, o++)
715 {
716 real_t kx = kv1[0]->GetBotella(ijk[0] + i);
717 if (!kv1[0]->inSpan(kx, ijk[0]+orders[0]+1)) { continue; }
718 ip.x = kv1[0]->GetRefPoint(kx, ijk[0]+orders[0]+1);
719
720 Trans.SetIntPoint(&ip);
721 vc.Eval(x, Trans, ip);
722
723 Trans.AdjugateJacobian().Mult(x,mx);
724 dofs(o) = mx(0);
725 }
726 }
727
728 for (int j = 0; j <= orders[1]+1; j++)
729 {
730 real_t ky = kv1[1]->GetBotella(ijk[1] + j);
731 if (!kv1[1]->inSpan(ky, ijk[1]+orders[1]+1))
732 {
733 o += orders[0] + 1;
734 continue;
735 }
736 ip.y = kv1[1]->GetRefPoint(ky, ijk[1]+orders[1]+1);
737 for (int i = 0; i <= orders[0]; i++, o++)
738 {
739 real_t kx = kv[0]->GetBotella(ijk[0] + i);
740 if (!kv[0]->inSpan(kx, ijk[0]+orders[0])) { continue; }
741 ip.x = kv[0]->GetRefPoint(kx, ijk[0]+orders[0]);
742
743 Trans.SetIntPoint(&ip);
744 vc.Eval(x, Trans, ip);
745
746 Trans.AdjugateJacobian().Mult(x,mx);
747 dofs(o) = mx(1);
748 }
749 }
750}
751
752NURBS_HDiv2DFiniteElement::~NURBS_HDiv2DFiniteElement()
753{
754 if (kv1[0]) { delete kv1[0]; }
755 if (kv1[1]) { delete kv1[1]; }
756}
757
758
759void NURBS_HDiv3DFiniteElement::SetOrder() const
760{
761 orders[0] = kv[0]->GetOrder();
762 orders[1] = kv[1]->GetOrder();
763 orders[2] = kv[2]->GetOrder();
764
765 if (kv1[0]) { delete kv1[0]; }
766 if (kv1[1]) { delete kv1[1]; }
767 if (kv1[2]) { delete kv1[2]; }
768
769 kv1[0] = kv[0]->DegreeElevate(1);
770 kv1[1] = kv[1]->DegreeElevate(1);
771 kv1[2] = kv[2]->DegreeElevate(1);
772
773 shape_x.SetSize(orders[0]+1);
774 shape_y.SetSize(orders[1]+1);
775 shape_z.SetSize(orders[2]+1);
776
777 dshape_x.SetSize(orders[0]+1);
778 dshape_y.SetSize(orders[1]+1);
779 dshape_z.SetSize(orders[2]+1);
780
781 d2shape_x.SetSize(orders[0]+1);
782 d2shape_y.SetSize(orders[1]+1);
783 d2shape_z.SetSize(orders[2]+1);
784
785 shape1_x.SetSize(orders[0]+2);
786 shape1_y.SetSize(orders[1]+2);
787 shape1_z.SetSize(orders[2]+2);
788
789 dshape1_x.SetSize(orders[0]+2);
790 dshape1_y.SetSize(orders[1]+2);
791 dshape1_z.SetSize(orders[2]+2);
792
793 d2shape1_x.SetSize(orders[0]+2);
794 d2shape1_y.SetSize(orders[1]+2);
795 d2shape1_z.SetSize(orders[2]+2);
796
797 order = max(orders[0]+1, max( orders[1]+1, orders[2]+1));
798 dof = (orders[0] + 2)*(orders[1] + 1)*(orders[2] + 1) +
799 (orders[0] + 1)*(orders[1] + 2)*(orders[2] + 1) +
800 (orders[0] + 1)*(orders[1] + 1)*(orders[2] + 2);
801 u.SetSize(dof);
802 du.SetSize(dof);
803 weights.SetSize(dof);
804}
805
806void NURBS_HDiv3DFiniteElement::CalcVShape(const IntegrationPoint &ip,
807 DenseMatrix &shape) const
808{
809 kv[0]->CalcShape(shape_x, ijk[0], ip.x);
810 kv[1]->CalcShape(shape_y, ijk[1], ip.y);
811 kv[2]->CalcShape(shape_z, ijk[2], ip.z);
812
813 kv1[0]->CalcShape(shape1_x, ijk[0], ip.x);
814 kv1[1]->CalcShape(shape1_y, ijk[1], ip.y);
815 kv1[2]->CalcShape(shape1_z, ijk[2], ip.z);
816
817 shape = 0.0;
818 int o = 0;
819 for (int k = 0; k <= orders[2]; k++)
820 {
821 const real_t sz = shape_z(k);
822 for (int j = 0; j <= orders[1]; j++)
823 {
824 const real_t sy_sz = shape_y(j)*sz;
825 for (int i = 0; i <= orders[0]+1; i++, o++)
826 {
827 shape(o,0) = shape1_x(i)*sy_sz;
828 }
829 }
830 }
831
832 for (int k = 0; k <= orders[2]; k++)
833 {
834 const real_t sz = shape_z(k);
835 for (int j = 0; j <= orders[1]+1; j++)
836 {
837 const real_t sy1_sz = shape1_y(j)*sz;
838 for (int i = 0; i <= orders[0]; i++, o++)
839 {
840 shape(o,1) = shape_x(i)*sy1_sz;
841 }
842 }
843 }
844
845 for (int k = 0; k <= orders[2]+1; k++)
846 {
847 const real_t sz1 = shape1_z(k);
848 for (int j = 0; j <= orders[1]; j++)
849 {
850 const real_t sy_sz1 = shape_y(j)*sz1;
851 for (int i = 0; i <= orders[0]; i++, o++)
852 {
853 shape(o,2) = shape_x(i)*sy_sz1;
854 }
855 }
856 }
857}
858
859void NURBS_HDiv3DFiniteElement::CalcVShape(ElementTransformation &Trans,
860 DenseMatrix &shape) const
861{
862 CalcVShape(Trans.GetIntPoint(), shape);
863 const DenseMatrix & J = Trans.Jacobian();
864 MFEM_ASSERT(J.Width() == 3 && J.Height() == 3,
865 "RT_R2D_FiniteElement cannot be embedded in "
866 "3 dimensional spaces");
867 for (int i=0; i<dof; i++)
868 {
869 real_t sx = shape(i, 0);
870 real_t sy = shape(i, 1);
871 real_t sz = shape(i, 2);
872 shape(i, 0) = sx * J(0, 0) + sy * J(0, 1) + sz * J(0, 2);
873 shape(i, 1) = sx * J(1, 0) + sy * J(1, 1) + sz * J(1, 2);
874 shape(i, 2) = sx * J(2, 0) + sy * J(2, 1) + sz * J(2, 2);
875 }
876 shape *= (1.0 / Trans.Weight());
877}
878
879void NURBS_HDiv3DFiniteElement::CalcDivShape(const IntegrationPoint &ip,
880 Vector &divshape) const
881{
882 kv[0]->CalcShape ( shape_x, ijk[0], ip.x);
883 kv[1]->CalcShape ( shape_y, ijk[1], ip.y);
884 kv[2]->CalcShape ( shape_z, ijk[2], ip.z);
885
886 kv1[0]->CalcDShape(dshape1_x, ijk[0], ip.x);
887 kv1[1]->CalcDShape(dshape1_y, ijk[1], ip.y);
888 kv1[2]->CalcDShape(dshape1_z, ijk[2], ip.z);
889
890 int o = 0;
891 for (int k = 0; k <= orders[2]; k++)
892 {
893 const real_t sz = shape_z(k);
894 for (int j = 0; j <= orders[1]; j++)
895 {
896 const real_t sy_sz = shape_y(j)*sz;
897 for (int i = 0; i <= orders[0]+1; i++, o++)
898 {
899 divshape(o) = dshape1_x(i)*sy_sz;
900 }
901 }
902 }
903
904 for (int k = 0; k <= orders[2]; k++)
905 {
906 const real_t sz = shape_z(k);
907 for (int j = 0; j <= orders[1]+1; j++)
908 {
909 const real_t dy1_sz = dshape1_y(j)*sz;
910 for (int i = 0; i <= orders[0]; i++, o++)
911 {
912 divshape(o) = shape_x(i)*dy1_sz;
913 }
914 }
915 }
916
917 for (int k = 0; k <= orders[2]+1; k++)
918 {
919 const real_t dz1 = dshape1_z(k);
920 for (int j = 0; j <= orders[1]; j++)
921 {
922 const real_t sy_dz1 = shape_y(j)*dz1;
923 for (int i = 0; i <= orders[0]; i++, o++)
924 {
925 divshape(o) = shape_x(i)*sy_dz1;
926 }
927 }
928 }
929}
930
931
932void NURBS_HDiv3DFiniteElement::Project(VectorCoefficient &vc,
933 ElementTransformation &Trans,
934 Vector &dofs) const
935{
936 MFEM_ASSERT(dofs.Size() == dof, "");
937 MFEM_ASSERT(vc.GetVDim() == 3, "");
938 Vector x(2), mx(3);
939 IntegrationPoint ip;
940
941 int o = 0;
942
943 for (int k = 0; k <= orders[2]; k++)
944 {
945 real_t kz = kv[2]->GetBotella(ijk[2] + k);
946 if (!kv[2]->inSpan(kz, ijk[2]+orders[2]))
947 {
948 o += (orders[0] + 2)*(orders[1] + 1);
949 continue;
950 }
951 ip.z = kv[2]->GetRefPoint(kz, ijk[2]+orders[2]);
952 for (int j = 0; j <= orders[1]; j++)
953 {
954 real_t ky = kv[1]->GetBotella(ijk[1] + j);
955 if (!kv[1]->inSpan(ky, ijk[1]+orders[1]))
956 {
957 o += orders[0] + 2;
958 continue;
959 }
960 ip.y = kv[1]->GetRefPoint(ky, ijk[1]+orders[1]);
961 for (int i = 0; i <= orders[0]+1; i++, o++)
962 {
963 real_t kx = kv1[0]->GetBotella(ijk[0] + i);
964 if (!kv1[0]->inSpan(kx, ijk[0]+orders[0]+1)) { continue; }
965 ip.x = kv1[0]->GetRefPoint(kx, ijk[0]+orders[0]+1);
966
967 Trans.SetIntPoint(&ip);
968 vc.Eval(x, Trans, ip);
969
970 Trans.AdjugateJacobian().Mult(x,mx);
971 dofs(o) = mx(0);
972 }
973 }
974 }
975
976 for (int k = 0; k <= orders[2]; k++)
977 {
978 real_t kz = kv[2]->GetBotella(ijk[2] + k);
979 if (!kv[2]->inSpan(kz, ijk[2]+orders[2]))
980 {
981 o += (orders[0] + 1)*(orders[1] + 2);
982 continue;
983 }
984 ip.z = kv[2]->GetRefPoint(kz, ijk[2]+orders[2]);
985 for (int j = 0; j <= orders[1]+1; j++)
986 {
987 real_t ky = kv1[1]->GetBotella(ijk[1] + j);
988 if (!kv1[1]->inSpan(ky, ijk[1]+orders[1]+1))
989 {
990 o += orders[0] + 1;
991 continue;
992 }
993 ip.y = kv1[1]->GetRefPoint(ky, ijk[1]+orders[1]+1);
994 for (int i = 0; i <= orders[0]; i++, o++)
995 {
996 real_t kx = kv[0]->GetBotella(ijk[0] + i);
997 if (!kv[0]->inSpan(kx, ijk[0]+orders[0])) { continue; }
998 ip.x = kv[0]->GetRefPoint(kx, ijk[0]+orders[0]);
999
1000 Trans.SetIntPoint(&ip);
1001 vc.Eval(x, Trans, ip);
1002
1003 Trans.AdjugateJacobian().Mult(x,mx);
1004 dofs(o) = mx(1);
1005 }
1006 }
1007 }
1008
1009 for (int k = 0; k <= orders[2]+1; k++)
1010 {
1011 real_t kz = kv1[2]->GetBotella(ijk[2] + k);
1012 if (!kv1[2]->inSpan(kz, ijk[2]+orders[2]+1))
1013 {
1014 o += (orders[0] + 1)*(orders[1] + 1);
1015 continue;
1016 }
1017 ip.z = kv1[2]->GetRefPoint(kz, ijk[2]+orders[2]+1);
1018 for (int j = 0; j <= orders[1]; j++)
1019 {
1020 real_t ky = kv[1]->GetBotella(ijk[1] + j);
1021 if (!kv[1]->inSpan(ky, ijk[1]+orders[1]))
1022 {
1023 o += orders[0] + 1;
1024 continue;
1025 }
1026 ip.y = kv[1]->GetRefPoint(ky, ijk[1]+orders[1]);
1027 for (int i = 0; i <= orders[0]; i++, o++)
1028 {
1029 real_t kx = kv[0]->GetBotella(ijk[0] + i);
1030 if (!kv[0]->inSpan(kx, ijk[0]+orders[0])) { continue; }
1031 ip.x = kv[0]->GetRefPoint(kx, ijk[0]+orders[0]);
1032
1033 Trans.SetIntPoint(&ip);
1034 vc.Eval(x, Trans, ip);
1035
1036 Trans.AdjugateJacobian().Mult(x,mx);
1037 dofs(o) = mx(2);
1038 }
1039 }
1040 }
1041
1042}
1043
1044
1045NURBS_HDiv3DFiniteElement::~NURBS_HDiv3DFiniteElement()
1046{
1047 if (kv1[0]) { delete kv1[0]; }
1048 if (kv1[1]) { delete kv1[1]; }
1049 if (kv1[2]) { delete kv1[2]; }
1050}
1051
1052void NURBS_HCurl2DFiniteElement::SetOrder() const
1053{
1054 orders[0] = kv[0]->GetOrder();
1055 orders[1] = kv[1]->GetOrder();
1056
1057 if (kv1[0]) { delete kv1[0]; }
1058 if (kv1[1]) { delete kv1[1]; }
1059
1060 kv1[0] = kv[0]->DegreeElevate(1);
1061 kv1[1] = kv[1]->DegreeElevate(1);
1062
1063 shape_x.SetSize(orders[0]+1);
1064 shape_y.SetSize(orders[1]+1);
1065
1066 dshape_x.SetSize(orders[0]+1);
1067 dshape_y.SetSize(orders[1]+1);
1068
1069 d2shape_x.SetSize(orders[0]+1);
1070 d2shape_y.SetSize(orders[1]+1);
1071
1072 shape1_x.SetSize(orders[0]+2);
1073 shape1_y.SetSize(orders[1]+2);
1074
1075 dshape1_x.SetSize(orders[0]+2);
1076 dshape1_y.SetSize(orders[1]+2);
1077
1078 d2shape1_x.SetSize(orders[0]+2);
1079 d2shape1_y.SetSize(orders[1]+2);
1080
1081 order = max(orders[0]+1, orders[1]+1);
1082 dof = (orders[0] + 1)*(orders[1] + 2)
1083 + (orders[1] + 2)*(orders[1] + 1);
1084 u.SetSize(dof);
1085 du.SetSize(dof);
1086 weights.SetSize(dof);
1087}
1088
1089void NURBS_HCurl2DFiniteElement::CalcVShape(const IntegrationPoint &ip,
1090 DenseMatrix &shape) const
1091{
1092 kv[0]->CalcShape(shape_x, ijk[0], ip.x);
1093 kv[1]->CalcShape(shape_y, ijk[1], ip.y);
1094
1095 kv1[0]->CalcShape(shape1_x, ijk[0], ip.x);
1096 kv1[1]->CalcShape(shape1_y, ijk[1], ip.y);
1097
1098 int o = 0;
1099 for (int j = 0; j <= orders[1]+1; j++)
1100 {
1101 const real_t sy1 = shape1_y(j);
1102 for (int i = 0; i <= orders[0]; i++, o++)
1103 {
1104 shape(o,0) = shape_x(i)*sy1;
1105 shape(o,1) = 0.0;
1106 }
1107 }
1108
1109 for (int j = 0; j <= orders[1]; j++)
1110 {
1111 const real_t sy = shape_y(j);
1112 for (int i = 0; i <= orders[0]+1; i++, o++)
1113 {
1114 shape(o,0) = 0.0;
1115 shape(o,1) = shape1_x(i)*sy;
1116 }
1117 }
1118}
1119
1120void NURBS_HCurl2DFiniteElement::CalcVShape(ElementTransformation &Trans,
1121 DenseMatrix &shape) const
1122{
1123 CalcVShape(Trans.GetIntPoint(), shape);
1124 const DenseMatrix & JI = Trans.InverseJacobian();
1125 MFEM_ASSERT(JI.Width() == 2 && JI.Height() == 2,
1126 "NURBS_HCurl2DFiniteElement cannot be embedded in "
1127 "3 dimensional spaces");
1128 for (int i=0; i<dof; i++)
1129 {
1130 real_t sx = shape(i, 0);
1131 real_t sy = shape(i, 1);
1132 shape(i, 0) = sx * JI(0, 0) + sy * JI(1, 0);
1133 shape(i, 1) = sx * JI(0, 1) + sy * JI(1, 1);
1134 }
1135}
1136
1137void NURBS_HCurl2DFiniteElement::CalcCurlShape(const IntegrationPoint &ip,
1138 DenseMatrix &curl_shape) const
1139{
1140 kv[0]->CalcShape ( shape_x, ijk[0], ip.x);
1141 kv[1]->CalcShape ( shape_y, ijk[1], ip.y);
1142
1143 kv1[0]->CalcDShape(dshape1_x, ijk[0], ip.x);
1144 kv1[1]->CalcDShape(dshape1_y, ijk[1], ip.y);
1145
1146 int o = 0;
1147 for (int j = 0; j <= orders[1]+1; j++)
1148 {
1149 const real_t dsy1 = dshape1_y(j);
1150 for (int i = 0; i <= orders[0]; i++, o++)
1151 {
1152 curl_shape(o,0) = -shape_x(i)*dsy1;
1153 }
1154 }
1155
1156 for (int j = 0; j <= orders[1]; j++)
1157 {
1158 const real_t sy = shape_y(j);
1159 for (int i = 0; i <= orders[0]+1; i++, o++)
1160 {
1161 curl_shape(o,0) = dshape1_x(i)*sy;
1162 }
1163 }
1164}
1165
1166void NURBS_HCurl2DFiniteElement::Project(VectorCoefficient &vc,
1167 ElementTransformation &Trans,
1168 Vector &dofs) const
1169{
1170 MFEM_ASSERT(dofs.Size() == dof, "");
1171 MFEM_ASSERT(vc.GetVDim() == 2, "");
1172 Vector x(2), xm(2);
1173 IntegrationPoint ip;
1174 int i, j, o;
1175 for (o = 0, j = 0; j <= orders[1]+1; j++)
1176 {
1177 real_t ky = kv1[1]->GetBotella(ijk[1] + j);
1178 if (!kv1[1]->inSpan(ky, ijk[1]+orders[1]+1))
1179 {
1180 o += orders[0] + 1;
1181 continue;
1182 }
1183 ip.y = kv1[1]->GetRefPoint(ky, ijk[1]+orders[1]+1);
1184 for (i = 0; i <= orders[0]; i++, o++)
1185 {
1186 real_t kx = kv[0]->GetBotella(ijk[0] + i);
1187 if (!kv[0]->inSpan(kx, ijk[0]+orders[0])) { continue; }
1188 ip.x = kv[0]->GetRefPoint(kx, ijk[0]+orders[0]);
1189
1190 Trans.SetIntPoint(&ip);
1191 vc.Eval(x, Trans, ip);
1192
1193 Trans.Jacobian().MultTranspose(x,xm);
1194 dofs(o) = xm(0);
1195 }
1196 }
1197
1198 for (j = 0; j <= orders[1]; j++)
1199 {
1200 real_t ky = kv[1]->GetBotella(ijk[1] + j);
1201 if (!kv[1]->inSpan(ky, ijk[1]+orders[1]))
1202 {
1203 o += orders[0] + 2;
1204 continue;
1205 }
1206 ip.y = kv[1]->GetRefPoint(ky, ijk[1]+orders[1]);
1207 for (i = 0; i <= orders[0]+1; i++, o++)
1208 {
1209 real_t kx = kv1[0]->GetBotella(ijk[0] + i);
1210 if (!kv1[0]->inSpan(kx, ijk[0]+orders[0]+1)) { continue; }
1211 ip.x = kv1[0]->GetRefPoint(kx, ijk[0]+orders[0]+1);
1212
1213 Trans.SetIntPoint(&ip);
1214 vc.Eval(x, Trans, ip);
1215
1216 Trans.Jacobian().MultTranspose(x,xm);
1217 dofs(o) = xm(1);
1218 }
1219 }
1220}
1221
1222NURBS_HCurl2DFiniteElement::~NURBS_HCurl2DFiniteElement()
1223{
1224 if (kv1[0]) { delete kv1[0]; }
1225 if (kv1[1]) { delete kv1[1]; }
1226}
1227
1228void NURBS_HCurl3DFiniteElement::SetOrder() const
1229{
1230 orders[0] = kv[0]->GetOrder();
1231 orders[1] = kv[1]->GetOrder();
1232 orders[2] = kv[2]->GetOrder();
1233
1234 if (kv1[0]) { delete kv1[0]; }
1235 if (kv1[1]) { delete kv1[1]; }
1236 if (kv1[2]) { delete kv1[2]; }
1237
1238 kv1[0] = kv[0]->DegreeElevate(1);
1239 kv1[1] = kv[1]->DegreeElevate(1);
1240 kv1[2] = kv[2]->DegreeElevate(1);
1241
1242 shape_x.SetSize(orders[0]+1);
1243 shape_y.SetSize(orders[1]+1);
1244 shape_z.SetSize(orders[2]+1);
1245
1246 dshape_x.SetSize(orders[0]+1);
1247 dshape_y.SetSize(orders[1]+1);
1248 dshape_z.SetSize(orders[2]+1);
1249
1250 d2shape_x.SetSize(orders[0]+1);
1251 d2shape_y.SetSize(orders[1]+1);
1252 d2shape_z.SetSize(orders[2]+1);
1253
1254 shape1_x.SetSize(orders[0]+2);
1255 shape1_y.SetSize(orders[1]+2);
1256 shape1_z.SetSize(orders[2]+2);
1257
1258 dshape1_x.SetSize(orders[0]+2);
1259 dshape1_y.SetSize(orders[1]+2);
1260 dshape1_z.SetSize(orders[2]+2);
1261
1262 d2shape1_x.SetSize(orders[0]+2);
1263 d2shape1_y.SetSize(orders[1]+2);
1264 d2shape1_z.SetSize(orders[2]+2);
1265
1266 order = max(orders[0]+1, max( orders[1]+1, orders[2]+1));
1267 dof = (orders[0] + 1)*(orders[1] + 2)*(orders[2] + 2) +
1268 (orders[0] + 2)*(orders[1] + 1)*(orders[2] + 2) +
1269 (orders[0] + 2)*(orders[1] + 2)*(orders[2] + 1);
1270 u.SetSize(dof);
1271 du.SetSize(dof);
1272 weights.SetSize(dof);
1273}
1274
1275void NURBS_HCurl3DFiniteElement::CalcVShape(const IntegrationPoint &ip,
1276 DenseMatrix &shape) const
1277{
1278 kv[0]->CalcShape(shape_x, ijk[0], ip.x);
1279 kv[1]->CalcShape(shape_y, ijk[1], ip.y);
1280 kv[2]->CalcShape(shape_z, ijk[2], ip.z);
1281
1282 kv1[0]->CalcShape(shape1_x, ijk[0], ip.x);
1283 kv1[1]->CalcShape(shape1_y, ijk[1], ip.y);
1284 kv1[2]->CalcShape(shape1_z, ijk[2], ip.z);
1285
1286 shape = 0.0;
1287 int o = 0;
1288 for (int k = 0; k <= orders[2]+1; k++)
1289 {
1290 const real_t sz1 = shape1_z(k);
1291 for (int j = 0; j <= orders[1]+1; j++)
1292 {
1293 const real_t sy1_sz1 = shape1_y(j)*sz1;
1294 for (int i = 0; i <= orders[0]; i++, o++)
1295 {
1296 shape(o,0) = shape_x(i)*sy1_sz1;
1297 }
1298 }
1299 }
1300
1301 for (int k = 0; k <= orders[2]+1; k++)
1302 {
1303 const real_t sz1 = shape1_z(k);
1304 for (int j = 0; j <= orders[1]; j++)
1305 {
1306 const real_t sy_sz1 = shape_y(j)*sz1;
1307 for (int i = 0; i <= orders[0]+1; i++, o++)
1308 {
1309 shape(o,1) = shape1_x(i)*sy_sz1;
1310 }
1311 }
1312 }
1313
1314 for (int k = 0; k <= orders[2]; k++)
1315 {
1316 const real_t sz = shape_z(k);
1317 for (int j = 0; j <= orders[1]+1; j++)
1318 {
1319 const real_t sy1_sz = shape1_y(j)*sz;
1320 for (int i = 0; i <= orders[0]+1; i++, o++)
1321 {
1322 shape(o,2) = shape1_x(i)*sy1_sz;
1323 }
1324 }
1325 }
1326}
1327
1328void NURBS_HCurl3DFiniteElement::CalcVShape(ElementTransformation &Trans,
1329 DenseMatrix &shape) const
1330{
1331 CalcVShape(Trans.GetIntPoint(), shape);
1332 const DenseMatrix & JI = Trans.InverseJacobian();
1333 MFEM_ASSERT(JI.Width() == 3 && JI.Height() == 3,
1334 "NURBS_HCurl3DFiniteElement must be in a"
1335 "3 dimensional spaces");
1336 for (int i=0; i<dof; i++)
1337 {
1338 real_t sx = shape(i, 0);
1339 real_t sy = shape(i, 1);
1340 real_t sz = shape(i, 2);
1341 shape(i, 0) = sx * JI(0, 0) + sy * JI(1, 0) + sz * JI(2, 0);
1342 shape(i, 1) = sx * JI(0, 1) + sy * JI(1, 1) + sz * JI(2, 1);
1343 shape(i, 2) = sx * JI(0, 2) + sy * JI(1, 2) + sz * JI(2, 2);
1344 }
1345}
1346
1347void NURBS_HCurl3DFiniteElement::CalcCurlShape(const IntegrationPoint &ip,
1348 DenseMatrix &curl_shape) const
1349{
1350 kv[0]->CalcShape ( shape_x, ijk[0], ip.x);
1351 kv[1]->CalcShape ( shape_y, ijk[1], ip.y);
1352 kv[2]->CalcShape ( shape_z, ijk[2], ip.z);
1353
1354 kv1[0]->CalcShape(shape1_x, ijk[0], ip.x);
1355 kv1[1]->CalcShape(shape1_y, ijk[1], ip.y);
1356 kv1[2]->CalcShape(shape1_z, ijk[2], ip.z);
1357
1358 kv1[0]->CalcDShape(dshape1_x, ijk[0], ip.x);
1359 kv1[1]->CalcDShape(dshape1_y, ijk[1], ip.y);
1360 kv1[2]->CalcDShape(dshape1_z, ijk[2], ip.z);
1361
1362 int o = 0;
1363 for (int k = 0; k <= orders[2]+1; k++)
1364 {
1365 const real_t sz1 = shape1_z(k), dsz1 = dshape1_z(k);
1366 for (int j = 0; j <= orders[1]+1; j++)
1367 {
1368 const real_t sy1_dsz1 = shape1_y(j)*dsz1,
1369 dsy1_sz1 = dshape1_y(j)*sz1;
1370 for (int i = 0; i <= orders[0]; i++, o++)
1371 {
1372 curl_shape(o,0) = 0.0;
1373 curl_shape(o,1) = shape_x(i)*sy1_dsz1;
1374 curl_shape(o,2) = -shape_x(i)*dsy1_sz1;
1375 }
1376 }
1377 }
1378
1379 for (int k = 0; k <= orders[2]+1; k++)
1380 {
1381 const real_t sz1 = shape1_z(k), dsz1 = dshape1_z(k);
1382 for (int j = 0; j <= orders[1]; j++)
1383 {
1384 const real_t sy_dsz1 = shape_y(j)*dsz1,
1385 sy_sz1 = shape_y(j)*sz1;
1386 for (int i = 0; i <= orders[0]+1; i++, o++)
1387 {
1388 curl_shape(o,0) = -shape1_x(i)*sy_dsz1;
1389 curl_shape(o,1) = 0.0;
1390 curl_shape(o,2) = dshape1_x(i)*sy_sz1;
1391 }
1392 }
1393 }
1394
1395 for (int k = 0; k <= orders[2]; k++)
1396 {
1397 const real_t sz = shape_z(k);
1398 for (int j = 0; j <= orders[1]+1; j++)
1399 {
1400 const real_t sy1_sz = shape1_y(j)*sz,
1401 dsy1_sz = dshape1_y(j)*sz;
1402 for (int i = 0; i <= orders[0]+1; i++, o++)
1403 {
1404 curl_shape(o,0) = shape1_x(i)*dsy1_sz;
1405 curl_shape(o,1) = -dshape1_x(i)*sy1_sz;
1406 curl_shape(o,2) = 0.0;
1407
1408 }
1409 }
1410 }
1411}
1412
1413void NURBS_HCurl3DFiniteElement::Project(VectorCoefficient &vc,
1414 ElementTransformation &Trans,
1415 Vector &dofs) const
1416{
1417 MFEM_ASSERT(dofs.Size() == dof, "");
1418 MFEM_ASSERT(vc.GetVDim() == 3, "");
1419 Vector x(3), xm(3);
1420 IntegrationPoint ip;
1421
1422 int o = 0;
1423 for (int k = 0; k <= orders[2]+1; k++)
1424 {
1425 real_t kz = kv1[2]->GetBotella(ijk[2] + k);
1426 if (!kv1[2]->inSpan(kz, ijk[2]+orders[2]+1))
1427 {
1428 o += (orders[0] + 1)*(orders[1] + 2);
1429 continue;
1430 }
1431 ip.z = kv1[2]->GetRefPoint(kz, ijk[2]+orders[2]+1);
1432 for (int j = 0; j <= orders[1]+1; j++)
1433 {
1434 real_t ky = kv1[1]->GetBotella(ijk[1] + j);
1435 if (!kv1[1]->inSpan(ky, ijk[1]+orders[1]+1))
1436 {
1437 o += orders[0] + 1;
1438 continue;
1439 }
1440 ip.y = kv1[1]->GetRefPoint(ky, ijk[1]+orders[1]+1);
1441 for (int i = 0; i <= orders[0]; i++, o++)
1442 {
1443 real_t kx = kv[0]->GetBotella(ijk[0] + i);
1444 if (!kv[0]->inSpan(kx, ijk[0]+orders[0])) { continue; }
1445 ip.x = kv[0]->GetRefPoint(kx, ijk[0]+orders[0]);
1446
1447 Trans.SetIntPoint(&ip);
1448 vc.Eval(x, Trans, ip);
1449
1450 Trans.Jacobian().MultTranspose(x,xm);
1451 dofs(o) = xm(0);
1452 }
1453 }
1454 }
1455
1456 for (int k = 0; k <= orders[2]+1; k++)
1457 {
1458 real_t kz = kv1[2]->GetBotella(ijk[2] + k);
1459 if (!kv1[2]->inSpan(kz, ijk[2]+orders[2]+1))
1460 {
1461 o += (orders[0] + 2)*(orders[1] + 1);
1462 continue;
1463 }
1464 ip.z = kv1[2]->GetRefPoint(kz, ijk[2]+orders[2]+1);
1465 for (int j = 0; j <= orders[1]; j++)
1466 {
1467 real_t ky = kv[1]->GetBotella(ijk[1] + j);
1468 if (!kv[1]->inSpan(ky, ijk[1]+orders[1]))
1469 {
1470 o += orders[0] + 2;
1471 continue;
1472 }
1473 ip.y = kv[1]->GetRefPoint(ky, ijk[1]+orders[1]);
1474 for (int i = 0; i <= orders[0]+1; i++, o++)
1475 {
1476 real_t kx = kv1[0]->GetBotella(ijk[0] + i);
1477 if (!kv1[0]->inSpan(kx, ijk[0]+orders[0]+1)) { continue; }
1478 ip.x = kv1[0]->GetRefPoint(kx, ijk[0]+orders[0]+1);
1479
1480 Trans.SetIntPoint(&ip);
1481 vc.Eval(x, Trans, ip);
1482
1483 Trans.Jacobian().MultTranspose(x,xm);
1484 dofs(o) = xm(1);
1485 }
1486 }
1487 }
1488
1489 for (int k = 0; k <= orders[2]; k++)
1490 {
1491 real_t kz = kv[2]->GetBotella(ijk[2] + k);
1492 if (!kv[2]->inSpan(kz, ijk[2]+orders[2]))
1493 {
1494 o += (orders[0] + 2)*(orders[1] + 2);
1495 continue;
1496 }
1497 ip.z = kv[2]->GetRefPoint(kz, ijk[2]+orders[2]);
1498 for (int j = 0; j <= orders[1]+1; j++)
1499 {
1500 real_t ky = kv1[1]->GetBotella(ijk[1] + j);
1501 if (!kv1[1]->inSpan(ky, ijk[1]+orders[1]+1))
1502 {
1503 o += orders[0] + 2;
1504 continue;
1505 }
1506 ip.y = kv1[1]->GetRefPoint(ky, ijk[1]+orders[1]+1);
1507 for (int i = 0; i <= orders[0]+1; i++, o++)
1508 {
1509 real_t kx = kv1[0]->GetBotella(ijk[0] + i);
1510 if (!kv1[0]->inSpan(kx, ijk[0]+orders[0]+1)) { continue; }
1511 ip.x = kv1[0]->GetRefPoint(kx, ijk[0]+orders[0]+1);
1512
1513 Trans.SetIntPoint(&ip);
1514 vc.Eval(x, Trans, ip);
1515
1516 Trans.Jacobian().MultTranspose(x,xm);
1517 dofs(o) = xm(2);
1518 }
1519 }
1520 }
1521
1522}
1523
1524
1525NURBS_HCurl3DFiniteElement::~NURBS_HCurl3DFiniteElement()
1526{
1527 if (kv1[0]) { delete kv1[0]; }
1528 if (kv1[1]) { delete kv1[1]; }
1529 if (kv1[2]) { delete kv1[2]; }
1530}
1531
1532}
mfem::Coefficient
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition coefficient.hpp:42
mfem::Coefficient::Eval
virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
mfem::DenseMatrix
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
mfem::DenseMatrix::Mult
void Mult(const real_t *x, real_t *y) const
Matrix vector multiplication.
Definition densemat.cpp:108
mfem::DenseMatrix::MultTranspose
void MultTranspose(const real_t *x, real_t *y) const
Multiply a vector with the transpose matrix.
Definition densemat.cpp:158
mfem::DenseMatrix::Data
real_t * Data() const
Returns the matrix data array. Warning: this method casts away constness.
Definition densemat.hpp:131
mfem::DenseMatrix::SetSize
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition densemat.hpp:125
mfem::ElementTransformation::InverseJacobian
const DenseMatrix & InverseJacobian()
Return the inverse of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition eltrans.hpp:158
mfem::ElementTransformation::AdjugateJacobian
const DenseMatrix & AdjugateJacobian()
Return the adjugate of the Jacobian matrix of the transformation at the currently set IntegrationPoin...
Definition eltrans.hpp:148
mfem::ElementTransformation::GetIntPoint
const IntegrationPoint & GetIntPoint()
Get a const reference to the currently set integration point. This will return NULL if no integration...
Definition eltrans.hpp:111
mfem::ElementTransformation::Weight
real_t Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint....
Definition eltrans.hpp:144
mfem::ElementTransformation::Jacobian
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint,...
Definition eltrans.hpp:132
mfem::ElementTransformation::SetIntPoint
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition eltrans.hpp:106
mfem::FiniteElement::dof
int dof
Number of degrees of freedom.
Definition fe_base.hpp:303
mfem::FiniteElement::orders
int orders[Geometry::MaxDim]
Anisotropic orders.
Definition fe_base.hpp:305
mfem::FiniteElement::order
int order
Order/degree of the shape functions.
Definition fe_base.hpp:304
mfem::IntegrationPoint
Class for integration point with weight.
Definition intrules.hpp:35
mfem::NURBS1DFiniteElement::CalcShape
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
mfem::NURBS1DFiniteElement::SetOrder
void SetOrder() const override
Definition fe_nurbs.cpp:22
mfem::NURBS1DFiniteElement::CalcHessian
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
mfem::NURBS1DFiniteElement::Project
void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const override
Definition fe_nurbs.cpp:87
mfem::NURBS1DFiniteElement::CalcDShape
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
mfem::NURBS2DFiniteElement::CalcDShape
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:165
mfem::NURBS2DFiniteElement::Project
void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const override
Definition fe_nurbs.cpp:258
mfem::NURBS2DFiniteElement::SetOrder
void SetOrder() const override
Definition fe_nurbs.cpp:128
mfem::NURBS2DFiniteElement::CalcHessian
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:200
mfem::NURBS2DFiniteElement::CalcShape
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:146
mfem::NURBS3DFiniteElement::SetOrder
void SetOrder() const override
Definition fe_nurbs.cpp:316
mfem::NURBS3DFiniteElement::CalcDShape
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:364
mfem::NURBS3DFiniteElement::CalcHessian
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:410
mfem::NURBS3DFiniteElement::Project
void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const override
Definition fe_nurbs.cpp:500
mfem::NURBS3DFiniteElement::CalcShape
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:340
mfem::NURBSFiniteElement::kv
Array< const KnotVector * > kv
Definition fe_nurbs.hpp:26
mfem::NURBS_HCurl2DFiniteElement::dshape1_y
Vector dshape1_y
Definition fe_nurbs.hpp:404
mfem::NURBS_HCurl2DFiniteElement::shape1_y
Vector shape1_y
Definition fe_nurbs.hpp:404
mfem::NURBS_HCurl2DFiniteElement::d2shape1_y
Vector d2shape1_y
Definition fe_nurbs.hpp:404
mfem::NURBS_HCurl2DFiniteElement::du
DenseMatrix du
Definition fe_nurbs.hpp:406
mfem::NURBS_HCurl2DFiniteElement::u
Vector u
Definition fe_nurbs.hpp:405
mfem::NURBS_HCurl2DFiniteElement::dshape1_x
Vector dshape1_x
Definition fe_nurbs.hpp:404
mfem::NURBS_HCurl2DFiniteElement::dshape_y
Vector dshape_y
Definition fe_nurbs.hpp:403
mfem::NURBS_HCurl2DFiniteElement::d2shape_y
Vector d2shape_y
Definition fe_nurbs.hpp:403
mfem::NURBS_HCurl2DFiniteElement::CalcCurlShape
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:1137
mfem::NURBS_HCurl2DFiniteElement::kv1
Array< const KnotVector * > kv1
Definition fe_nurbs.hpp:407
mfem::NURBS_HCurl2DFiniteElement::~NURBS_HCurl2DFiniteElement
~NURBS_HCurl2DFiniteElement()
Definition fe_nurbs.cpp:1222
mfem::NURBS_HCurl2DFiniteElement::d2shape_x
Vector d2shape_x
Definition fe_nurbs.hpp:403
mfem::NURBS_HCurl2DFiniteElement::shape_y
Vector shape_y
Definition fe_nurbs.hpp:403
mfem::NURBS_HCurl2DFiniteElement::Project
void Project(VectorCoefficient &vcoeff, ElementTransformation &Trans, Vector &dofs) const override
Definition fe_nurbs.cpp:1166
mfem::NURBS_HCurl2DFiniteElement::dshape_x
Vector dshape_x
Definition fe_nurbs.hpp:403
mfem::NURBS_HCurl2DFiniteElement::CalcVShape
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:1089
mfem::NURBS_HCurl2DFiniteElement::SetOrder
void SetOrder() const override
Definition fe_nurbs.cpp:1052
mfem::NURBS_HCurl2DFiniteElement::shape1_x
Vector shape1_x
Definition fe_nurbs.hpp:404
mfem::NURBS_HCurl2DFiniteElement::d2shape1_x
Vector d2shape1_x
Definition fe_nurbs.hpp:404
mfem::NURBS_HCurl2DFiniteElement::shape_x
Vector shape_x
Definition fe_nurbs.hpp:403
mfem::NURBS_HCurl3DFiniteElement::shape_x
Vector shape_x
Definition fe_nurbs.hpp:489
mfem::NURBS_HCurl3DFiniteElement::d2shape_y
Vector d2shape_y
Definition fe_nurbs.hpp:491
mfem::NURBS_HCurl3DFiniteElement::shape_y
Vector shape_y
Definition fe_nurbs.hpp:489
mfem::NURBS_HCurl3DFiniteElement::dshape1_y
Vector dshape1_y
Definition fe_nurbs.hpp:493
mfem::NURBS_HCurl3DFiniteElement::shape1_y
Vector shape1_y
Definition fe_nurbs.hpp:492
mfem::NURBS_HCurl3DFiniteElement::d2shape1_y
Vector d2shape1_y
Definition fe_nurbs.hpp:494
mfem::NURBS_HCurl3DFiniteElement::~NURBS_HCurl3DFiniteElement
~NURBS_HCurl3DFiniteElement()
Definition fe_nurbs.cpp:1525
mfem::NURBS_HCurl3DFiniteElement::dshape_y
Vector dshape_y
Definition fe_nurbs.hpp:490
mfem::NURBS_HCurl3DFiniteElement::CalcVShape
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:1275
mfem::NURBS_HCurl3DFiniteElement::shape1_z
Vector shape1_z
Definition fe_nurbs.hpp:492
mfem::NURBS_HCurl3DFiniteElement::u
Vector u
Definition fe_nurbs.hpp:495
mfem::NURBS_HCurl3DFiniteElement::d2shape_z
Vector d2shape_z
Definition fe_nurbs.hpp:491
mfem::NURBS_HCurl3DFiniteElement::Project
void Project(VectorCoefficient &vcoeff, ElementTransformation &Trans, Vector &dofs) const override
Definition fe_nurbs.cpp:1413
mfem::NURBS_HCurl3DFiniteElement::dshape1_z
Vector dshape1_z
Definition fe_nurbs.hpp:493
mfem::NURBS_HCurl3DFiniteElement::shape_z
Vector shape_z
Definition fe_nurbs.hpp:489
mfem::NURBS_HCurl3DFiniteElement::SetOrder
void SetOrder() const override
Definition fe_nurbs.cpp:1228
mfem::NURBS_HCurl3DFiniteElement::d2shape1_z
Vector d2shape1_z
Definition fe_nurbs.hpp:494
mfem::NURBS_HCurl3DFiniteElement::dshape1_x
Vector dshape1_x
Definition fe_nurbs.hpp:493
mfem::NURBS_HCurl3DFiniteElement::d2shape1_x
Vector d2shape1_x
Definition fe_nurbs.hpp:494
mfem::NURBS_HCurl3DFiniteElement::du
DenseMatrix du
Definition fe_nurbs.hpp:496
mfem::NURBS_HCurl3DFiniteElement::shape1_x
Vector shape1_x
Definition fe_nurbs.hpp:492
mfem::NURBS_HCurl3DFiniteElement::dshape_x
Vector dshape_x
Definition fe_nurbs.hpp:490
mfem::NURBS_HCurl3DFiniteElement::CalcCurlShape
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:1347
mfem::NURBS_HCurl3DFiniteElement::dshape_z
Vector dshape_z
Definition fe_nurbs.hpp:490
mfem::NURBS_HCurl3DFiniteElement::kv1
Array< const KnotVector * > kv1
Definition fe_nurbs.hpp:497
mfem::NURBS_HCurl3DFiniteElement::d2shape_x
Vector d2shape_x
Definition fe_nurbs.hpp:491
mfem::NURBS_HDiv2DFiniteElement::u
Vector u
Definition fe_nurbs.hpp:220
mfem::NURBS_HDiv2DFiniteElement::CalcVShape
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:617
mfem::NURBS_HDiv2DFiniteElement::dshape1_y
Vector dshape1_y
Definition fe_nurbs.hpp:219
mfem::NURBS_HDiv2DFiniteElement::du
DenseMatrix du
Definition fe_nurbs.hpp:221
mfem::NURBS_HDiv2DFiniteElement::dshape_y
Vector dshape_y
Definition fe_nurbs.hpp:218
mfem::NURBS_HDiv2DFiniteElement::~NURBS_HDiv2DFiniteElement
~NURBS_HDiv2DFiniteElement()
Definition fe_nurbs.cpp:752
mfem::NURBS_HDiv2DFiniteElement::shape1_y
Vector shape1_y
Definition fe_nurbs.hpp:219
mfem::NURBS_HDiv2DFiniteElement::CalcDivShape
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:666
mfem::NURBS_HDiv2DFiniteElement::SetOrder
void SetOrder() const override
Definition fe_nurbs.cpp:580
mfem::NURBS_HDiv2DFiniteElement::shape_x
Vector shape_x
Definition fe_nurbs.hpp:218
mfem::NURBS_HDiv2DFiniteElement::dshape1_x
Vector dshape1_x
Definition fe_nurbs.hpp:219
mfem::NURBS_HDiv2DFiniteElement::d2shape1_x
Vector d2shape1_x
Definition fe_nurbs.hpp:219
mfem::NURBS_HDiv2DFiniteElement::shape1_x
Vector shape1_x
Definition fe_nurbs.hpp:219
mfem::NURBS_HDiv2DFiniteElement::d2shape_x
Vector d2shape_x
Definition fe_nurbs.hpp:218
mfem::NURBS_HDiv2DFiniteElement::d2shape_y
Vector d2shape_y
Definition fe_nurbs.hpp:218
mfem::NURBS_HDiv2DFiniteElement::kv1
Array< const KnotVector * > kv1
Definition fe_nurbs.hpp:222
mfem::NURBS_HDiv2DFiniteElement::dshape_x
Vector dshape_x
Definition fe_nurbs.hpp:218
mfem::NURBS_HDiv2DFiniteElement::shape_y
Vector shape_y
Definition fe_nurbs.hpp:218
mfem::NURBS_HDiv2DFiniteElement::d2shape1_y
Vector d2shape1_y
Definition fe_nurbs.hpp:219
mfem::NURBS_HDiv2DFiniteElement::Project
void Project(VectorCoefficient &vcoeff, ElementTransformation &Trans, Vector &dofs) const override
Definition fe_nurbs.cpp:695
mfem::NURBS_HDiv3DFiniteElement::d2shape1_z
Vector d2shape1_z
Definition fe_nurbs.hpp:312
mfem::NURBS_HDiv3DFiniteElement::~NURBS_HDiv3DFiniteElement
~NURBS_HDiv3DFiniteElement()
Definition fe_nurbs.cpp:1045
mfem::NURBS_HDiv3DFiniteElement::dshape1_x
Vector dshape1_x
Definition fe_nurbs.hpp:311
mfem::NURBS_HDiv3DFiniteElement::SetOrder
void SetOrder() const override
Definition fe_nurbs.cpp:759
mfem::NURBS_HDiv3DFiniteElement::u
Vector u
Definition fe_nurbs.hpp:313
mfem::NURBS_HDiv3DFiniteElement::d2shape1_y
Vector d2shape1_y
Definition fe_nurbs.hpp:312
mfem::NURBS_HDiv3DFiniteElement::shape_y
Vector shape_y
Definition fe_nurbs.hpp:307
mfem::NURBS_HDiv3DFiniteElement::shape1_x
Vector shape1_x
Definition fe_nurbs.hpp:310
mfem::NURBS_HDiv3DFiniteElement::d2shape_x
Vector d2shape_x
Definition fe_nurbs.hpp:309
mfem::NURBS_HDiv3DFiniteElement::dshape1_z
Vector dshape1_z
Definition fe_nurbs.hpp:311
mfem::NURBS_HDiv3DFiniteElement::kv1
Array< const KnotVector * > kv1
Definition fe_nurbs.hpp:315
mfem::NURBS_HDiv3DFiniteElement::dshape_x
Vector dshape_x
Definition fe_nurbs.hpp:308
mfem::NURBS_HDiv3DFiniteElement::dshape_y
Vector dshape_y
Definition fe_nurbs.hpp:308
mfem::NURBS_HDiv3DFiniteElement::CalcDivShape
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:879
mfem::NURBS_HDiv3DFiniteElement::shape1_y
Vector shape1_y
Definition fe_nurbs.hpp:310
mfem::NURBS_HDiv3DFiniteElement::d2shape_z
Vector d2shape_z
Definition fe_nurbs.hpp:309
mfem::NURBS_HDiv3DFiniteElement::shape1_z
Vector shape1_z
Definition fe_nurbs.hpp:310
mfem::NURBS_HDiv3DFiniteElement::dshape1_y
Vector dshape1_y
Definition fe_nurbs.hpp:311
mfem::NURBS_HDiv3DFiniteElement::du
DenseMatrix du
Definition fe_nurbs.hpp:314
mfem::NURBS_HDiv3DFiniteElement::Project
void Project(VectorCoefficient &vcoeff, ElementTransformation &Trans, Vector &dofs) const override
Definition fe_nurbs.cpp:932
mfem::NURBS_HDiv3DFiniteElement::shape_z
Vector shape_z
Definition fe_nurbs.hpp:307
mfem::NURBS_HDiv3DFiniteElement::dshape_z
Vector dshape_z
Definition fe_nurbs.hpp:308
mfem::NURBS_HDiv3DFiniteElement::d2shape_y
Vector d2shape_y
Definition fe_nurbs.hpp:309
mfem::NURBS_HDiv3DFiniteElement::shape_x
Vector shape_x
Definition fe_nurbs.hpp:307
mfem::NURBS_HDiv3DFiniteElement::CalcVShape
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:806
mfem::NURBS_HDiv3DFiniteElement::d2shape1_x
Vector d2shape1_x
Definition fe_nurbs.hpp:312
mfem::Operator::Height
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition operator.hpp:68
mfem::Operator::Width
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition operator.hpp:74
mfem::VectorCoefficient
Base class for vector Coefficients that optionally depend on time and space.
Definition coefficient.hpp:577
mfem::VectorCoefficient::GetVDim
int GetVDim()
Returns dimension of the vector.
Definition coefficient.hpp:593
mfem::VectorCoefficient::Eval
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...
mfem::Vector
Vector data type.
Definition vector.hpp:82
mfem::Vector::Size
int Size() const
Returns the size of the vector.
Definition vector.hpp:234
mfem::Vector::SetSize
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:633
fe_nurbs.hpp
mfem::add
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition vector.cpp:414
mfem::real_t
float real_t
Definition config.hpp:46
std
STL namespace.
Generated on Tue Sep 1 2026 16:39:57 for MFEM by doxygen 1.11.0