34 kv[0]->CalcShape(shape,
ijk[0], ip.
x);
37 for (
int i = 0; i <=
order; i++)
39 sum += (shape(i) *=
weights(i));
51 kv[0]->CalcDShape(grad,
ijk[0], ip.
x);
53 real_t sum = 0.0, dsum = 0.0;
54 for (
int i = 0; i <=
order; i++)
57 dsum += ( grad(i) *=
weights(i));
71 kv[0]->CalcDShape(grad,
ijk[0], ip.
x);
72 kv[0]->CalcD2Shape(hess,
ijk[0], ip.
x);
74 real_t sum = 0.0, dsum = 0.0, d2sum = 0.0;
75 for (
int i = 0; i <=
order; i++)
78 dsum += ( grad(i) *=
weights(i));
79 d2sum += ( hess(i) *=
weights(i));
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);
113 for (
int o = 0, j = 0; j <=
orders[1]; j++)
116 for (
int i = 0; i <=
orders[0]; i++, o++)
136 sum = dsum[0] = dsum[1] = 0.0;
137 for (
int o = 0, j = 0; j <=
orders[1]; j++)
140 for (
int i = 0; i <=
orders[0]; i++, o++)
153 for (
int o = 0; o <
dof; o++)
155 dshape(o,0) = dshape(o,0)*sum -
u(o)*dsum[0];
156 dshape(o,1) = dshape(o,1)*sum -
u(o)*dsum[1];
163 real_t sum, dsum[2], d2sum[3];
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++)
179 for (
int i = 0; i <=
orders[0]; i++, o++)
184 dsum[0] += (
du(o,0) = dsx*sy*
weights(o) );
185 dsum[1] += (
du(o,1) = sx*dsy*
weights(o) );
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) );
201 for (
int o = 0; o <
dof; o++)
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]);
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]);
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]);
251 for (
int o = 0, k = 0; k <=
orders[2]; k++)
254 for (
int j = 0; j <=
orders[1]; j++)
257 for (
int i = 0; i <=
orders[0]; i++, o++)
280 sum = dsum[0] = dsum[1] = dsum[2] = 0.0;
281 for (
int o = 0, k = 0; k <=
orders[2]; k++)
284 for (
int j = 0; j <=
orders[1]; j++)
289 for (
int i = 0; i <=
orders[0]; i++, o++)
305 for (
int o = 0; o <
dof; o++)
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];
316 real_t sum, dsum[3], d2sum[6];
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;
333 for (
int o = 0, k = 0; k <=
orders[2]; k++)
336 for (
int j = 0; j <=
orders[1]; j++)
339 for (
int i = 0; i <=
orders[0]; i++, o++)
342 sum += (
u(o) = sx*sy*sz*
weights(o) );
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) );
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) );
352 d2sum[3] += ( hessian(o,3) = sx*dsy*dsz*
weights(o) );
354 d2sum[4] += ( hessian(o,4) = sx*sy*d2sz*
weights(o) );
355 d2sum[5] += ( hessian(o,5) = sx*d2sy*sz*
weights(o) );
373 for (
int o = 0; o <
dof; o++)
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]);
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]);
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]);
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]);
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]);
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]);
410 if (
kv1[0]) {
delete kv1[0]; }
411 if (
kv1[1]) {
delete kv1[1]; }
413 kv1[0] =
kv[0]->DegreeElevate(1);
414 kv1[1] =
kv[1]->DegreeElevate(1);
452 for (
int j = 0; j <=
orders[1]; j++)
455 for (
int i = 0; i <=
orders[0]+1; i++, o++)
462 for (
int j = 0; j <=
orders[1]+1; j++)
465 for (
int i = 0; i <=
orders[0]; i++, o++)
479 "NURBS_HDiv2DFiniteElement cannot be embedded in "
480 "3 dimensional spaces");
481 for (
int i=0; i<
dof; i++)
485 shape(i, 0) = sx * J(0, 0) + sy * J(0, 1);
486 shape(i, 1) = sx * J(1, 0) + sy * J(1, 1);
488 shape *= (1.0 / Trans.
Weight());
501 for (
int j = 0; j <=
orders[1]; j++)
504 for (
int i = 0; i <=
orders[0]+1; i++, o++)
510 for (
int j = 0; j <=
orders[1]+1; j++)
513 for (
int i = 0; i <=
orders[0]; i++, o++)
522 if (
kv1[0]) {
delete kv1[0]; }
523 if (
kv1[1]) {
delete kv1[1]; }
533 if (
kv1[0]) {
delete kv1[0]; }
534 if (
kv1[1]) {
delete kv1[1]; }
535 if (
kv1[2]) {
delete kv1[2]; }
537 kv1[0] =
kv[0]->DegreeElevate(1);
538 kv1[1] =
kv[1]->DegreeElevate(1);
539 kv1[2] =
kv[2]->DegreeElevate(1);
587 for (
int k = 0; k <=
orders[2]; k++)
590 for (
int j = 0; j <=
orders[1]; j++)
593 for (
int i = 0; i <=
orders[0]+1; i++, o++)
600 for (
int k = 0; k <=
orders[2]; k++)
603 for (
int j = 0; j <=
orders[1]+1; j++)
606 for (
int i = 0; i <=
orders[0]; i++, o++)
608 shape(o,1) =
shape_x(i)*sy1_sz;
613 for (
int k = 0; k <=
orders[2]+1; k++)
616 for (
int j = 0; j <=
orders[1]; j++)
619 for (
int i = 0; i <=
orders[0]; i++, o++)
621 shape(o,2) =
shape_x(i)*sy_sz1;
633 "RT_R2D_FiniteElement cannot be embedded in "
634 "3 dimensional spaces");
635 for (
int i=0; i<
dof; i++)
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);
644 shape *= (1.0 / Trans.
Weight());
659 for (
int k = 0; k <=
orders[2]; k++)
662 for (
int j = 0; j <=
orders[1]; j++)
665 for (
int i = 0; i <=
orders[0]+1; i++, o++)
672 for (
int k = 0; k <=
orders[2]; k++)
675 for (
int j = 0; j <=
orders[1]+1; j++)
678 for (
int i = 0; i <=
orders[0]; i++, o++)
680 divshape(o) =
shape_x(i)*dy1_sz;
685 for (
int k = 0; k <=
orders[2]+1; k++)
688 for (
int j = 0; j <=
orders[1]; j++)
691 for (
int i = 0; i <=
orders[0]; i++, o++)
693 divshape(o) =
shape_x(i)*sy_dz1;
701 if (
kv1[0]) {
delete kv1[0]; }
702 if (
kv1[1]) {
delete kv1[1]; }
703 if (
kv1[2]) {
delete kv1[2]; }
711 if (
kv1[0]) {
delete kv1[0]; }
712 if (
kv1[1]) {
delete kv1[1]; }
714 kv1[0] =
kv[0]->DegreeElevate(1);
715 kv1[1] =
kv[1]->DegreeElevate(1);
753 for (
int j = 0; j <=
orders[1]+1; j++)
756 for (
int i = 0; i <=
orders[0]; i++, o++)
763 for (
int j = 0; j <=
orders[1]; j++)
766 for (
int i = 0; i <=
orders[0]+1; i++, o++)
780 "NURBS_HCurl2DFiniteElement cannot be embedded in "
781 "3 dimensional spaces");
782 for (
int i=0; i<
dof; i++)
786 shape(i, 0) = sx * JI(0, 0) + sy * JI(1, 0);
787 shape(i, 1) = sx * JI(0, 1) + sy * JI(1, 1);
801 for (
int j = 0; j <=
orders[1]+1; j++)
804 for (
int i = 0; i <=
orders[0]; i++, o++)
806 curl_shape(o,0) = -
shape_x(i)*dsy1;
810 for (
int j = 0; j <=
orders[1]; j++)
813 for (
int i = 0; i <=
orders[0]+1; i++, o++)
822 if (
kv1[0]) {
delete kv1[0]; }
823 if (
kv1[1]) {
delete kv1[1]; }
833 if (
kv1[0]) {
delete kv1[0]; }
834 if (
kv1[1]) {
delete kv1[1]; }
835 if (
kv1[2]) {
delete kv1[2]; }
837 kv1[0] =
kv[0]->DegreeElevate(1);
838 kv1[1] =
kv[1]->DegreeElevate(1);
839 kv1[2] =
kv[2]->DegreeElevate(1);
887 for (
int k = 0; k <=
orders[2]+1; k++)
890 for (
int j = 0; j <=
orders[1]+1; j++)
893 for (
int i = 0; i <=
orders[0]; i++, o++)
895 shape(o,0) =
shape_x(i)*sy1_sz1;
900 for (
int k = 0; k <=
orders[2]+1; k++)
903 for (
int j = 0; j <=
orders[1]; j++)
906 for (
int i = 0; i <=
orders[0]+1; i++, o++)
913 for (
int k = 0; k <=
orders[2]; k++)
916 for (
int j = 0; j <=
orders[1]+1; j++)
919 for (
int i = 0; i <=
orders[0]+1; i++, o++)
933 "NURBS_HCurl3DFiniteElement must be in a"
934 "3 dimensional spaces");
935 for (
int i=0; i<
dof; i++)
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);
962 for (
int k = 0; k <=
orders[2]+1; k++)
965 for (
int j = 0; j <=
orders[1]+1; j++)
969 for (
int i = 0; i <=
orders[0]; i++, o++)
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;
978 for (
int k = 0; k <=
orders[2]+1; k++)
981 for (
int j = 0; j <=
orders[1]; j++)
985 for (
int i = 0; i <=
orders[0]+1; i++, o++)
987 curl_shape(o,0) = -
shape1_x(i)*sy_dsz1;
988 curl_shape(o,1) = 0.0;
994 for (
int k = 0; k <=
orders[2]; k++)
997 for (
int j = 0; j <=
orders[1]+1; j++)
1001 for (
int i = 0; i <=
orders[0]+1; i++, o++)
1003 curl_shape(o,0) =
shape1_x(i)*dsy1_sz;
1005 curl_shape(o,2) = 0.0;
1013 if (
kv1[0]) {
delete kv1[0]; }
1014 if (
kv1[1]) {
delete kv1[1]; }
1015 if (
kv1[2]) {
delete kv1[2]; }
Data type dense matrix using column-major storage.
real_t * Data() const
Returns the matrix data array.
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
int dof
Number of degrees of freedom.
int orders[Geometry::MaxDim]
Anisotropic orders.
int order
Order/degree of the shape functions.
Class for integration point with weight.
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...
void SetOrder() const override
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...
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...
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...
void SetOrder() const override
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...
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...
void SetOrder() const override
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...
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...
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...
Array< const KnotVector * > kv
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...
Array< const KnotVector * > kv1
~NURBS_HCurl2DFiniteElement()
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...
void SetOrder() const override
~NURBS_HCurl3DFiniteElement()
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...
void SetOrder() const override
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...
Array< const KnotVector * > kv1
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...
~NURBS_HDiv2DFiniteElement()
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...
void SetOrder() const override
Array< const KnotVector * > kv1
~NURBS_HDiv3DFiniteElement()
void SetOrder() const override
Array< const KnotVector * > kv1
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...
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...
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
void SetSize(int s)
Resize the vector to size s.
void add(const Vector &v1, const Vector &v2, Vector &v)