16 void calcMapFromIdealElement(
int type, T &dSVec_dX, T &dSVec_dY, T &dSVec_dZ)
28 static const double cTri[2] = {-1. / std::sqrt(3.), 2. / std::sqrt(3.)};
29 dSVec_dY.scale(
cTri[1]);
30 dSVec_dY.axpy(dSVec_dX,
cTri[0]);
43 static const double cPyr = sqrt(2.);
49 static const double cTet[3] = {-3. / 2 / std::sqrt(6.),
50 -1. / 2 / std::sqrt(2.), std::sqrt(1.5)};
51 dSVec_dZ.scale(
cTet[2]);
52 dSVec_dZ.axpy(dSVec_dX,
cTet[0]);
53 dSVec_dZ.axpy(dSVec_dY,
cTet[1]);
60 inline double calcDet3D(
double M11,
double M12,
double M13,
double M21,
61 double M22,
double M23,
double M31,
double M32,
64 return M11 * (M22 * M33 - M23 * M32) - M12 * (M21 * M33 - M23 * M31) +
65 M13 * (M21 * M32 - M22 * M31);
71 inline void calcJDJ1D(
double dxdX,
double dxdY,
double dxdZ,
double dydX,
72 double dydY,
double dydZ,
double dzdX,
double dzdY,
73 double dzdZ,
int i,
int numMapNodes,
77 for(
int j = 0; j < numMapNodes; j++) {
78 const double &dPhidX = dShapeMat_dX(i, j);
79 JDJ(i, j) = dPhidX * (dydY * dzdZ - dzdY * dydZ);
80 JDJ(i, j + numMapNodes) = dPhidX * (dzdY * dxdZ - dxdY * dzdZ);
81 JDJ(i, j + 2 * numMapNodes) = dPhidX * (dxdY * dydZ - dydY * dxdZ);
83 JDJ(i, 3 * numMapNodes) =
84 calcDet3D(dxdX, dxdY, dxdZ, dydX, dydY, dydZ, dzdX, dzdY, dzdZ);
90 inline void calcJDJ2D(
double dxdX,
double dxdY,
double dxdZ,
double dydX,
91 double dydY,
double dydZ,
double dzdX,
double dzdY,
92 double dzdZ,
int i,
int numMapNodes,
97 for(
int j = 0; j < numMapNodes; j++) {
98 const double &dPhidX = dShapeMat_dX(i, j);
99 const double &dPhidY = dShapeMat_dY(i, j);
100 JDJ(i, j) = dPhidX * (dydY * dzdZ - dzdY * dydZ) +
101 dPhidY * (dzdX * dydZ - dydX * dzdZ);
102 JDJ(i, j + numMapNodes) = dPhidX * (dzdY * dxdZ - dxdY * dzdZ) +
103 dPhidY * (dxdX * dzdZ - dzdX * dxdZ);
104 JDJ(i, j + 2 * numMapNodes) = dPhidX * (dxdY * dydZ - dydY * dxdZ) +
105 dPhidY * (dydX * dxdZ - dxdX * dydZ);
107 JDJ(i, 3 * numMapNodes) =
108 calcDet3D(dxdX, dxdY, dxdZ, dydX, dydY, dydZ, dzdX, dzdY, dzdZ);
113 inline void calcJDJ3D(
double dxdX,
double dxdY,
double dxdZ,
double dydX,
114 double dydY,
double dydZ,
double dzdX,
double dzdY,
115 double dzdZ,
int i,
int numMapNodes,
121 for(
int j = 0; j < numMapNodes; j++) {
122 const double &dPhidX = dShapeMat_dX(i, j);
123 const double &dPhidY = dShapeMat_dY(i, j);
124 const double &dPhidZ = dShapeMat_dZ(i, j);
125 JDJ(i, j) = dPhidX * (dydY * dzdZ - dzdY * dydZ) +
126 dPhidY * (dzdX * dydZ - dydX * dzdZ) +
127 dPhidZ * (dydX * dzdY - dzdX * dydY);
128 JDJ(i, j + numMapNodes) = dPhidX * (dzdY * dxdZ - dxdY * dzdZ) +
129 dPhidY * (dxdX * dzdZ - dzdX * dxdZ) +
130 dPhidZ * (dzdX * dxdY - dxdX * dzdY);
131 JDJ(i, j + 2 * numMapNodes) = dPhidX * (dxdY * dydZ - dydY * dxdZ) +
132 dPhidY * (dydX * dxdZ - dxdX * dydZ) +
133 dPhidZ * (dxdX * dydY - dydX * dxdY);
135 JDJ(i, 3 * numMapNodes) =
136 calcDet3D(dxdX, dxdY, dxdZ, dydX, dydY, dydZ, dzdX, dzdY, dzdZ);
147 : _elementTag(elementTag), _data(fsdata)
158 const int numSampPnts = samplingPoints.
size1();
163 mapBasis->
df(samplingPoints, allDPsi);
164 const int numMapNodes = allDPsi.
size2();
169 for(
int i = 0; i < numSampPnts; i++) {
170 for(
int j = 0; j < numMapNodes; j++) {
235 calcMapFromIdealElement(type, dSMat_dX, dSMat_dY, dSMat_dZ);
242 calcMapFromIdealElement(type, dSVec_dX, dSVec_dY, dSVec_dZ);
248 dxyzdZ(jac[2], 1, 3);
253 : _elementTag(elementTag), _data(data), _dim(data.
getDimension())
255 const int parentType = data.
getType();
279 double xBar = 0., yBar = 0., zBar = 0.;
280 double barycenter[3] = {0., 0., 0.};
282 for(
int j = 0; j < primMapBasis->
points.
size2(); ++j) {
283 barycenter[j] += primMapBasis->
points(i, j);
291 primMapBasis->
df(xBar, yBar, zBar, barDPsi);
322 mapBasis->
df(lagPointsFast, allDPsiFast);
349 if((fabs(dxyzdXbar(0)) >= fabs(dxyzdXbar(1)) &&
350 fabs(dxyzdXbar(0)) >= fabs(dxyzdXbar(2))) ||
351 (fabs(dxyzdXbar(1)) >= fabs(dxyzdXbar(0)) &&
352 fabs(dxyzdXbar(1)) >= fabs(dxyzdXbar(2)))) {
353 result(0, 0) = dxyzdXbar(1);
354 result(0, 1) = -dxyzdXbar(0);
359 result(0, 1) = dxyzdXbar(2);
360 result(0, 2) = -dxyzdXbar(1);
363 sqrt(result(0, 0) * result(0, 0) + result(0, 1) * result(0, 1) +
364 result(0, 2) * result(0, 2));
365 result(0, 0) /= norm0;
366 result(0, 1) /= norm0;
367 result(0, 2) /= norm0;
369 result(1, 2) = dxyzdXbar(0) * result(0, 1) - dxyzdXbar(1) * result(0, 0);
370 result(1, 1) = -dxyzdXbar(0) * result(0, 2) + dxyzdXbar(2) * result(0, 0);
371 result(1, 0) = dxyzdXbar(1) * result(0, 2) - dxyzdXbar(2) * result(0, 1);
373 sqrt(result(1, 0) * result(1, 0) + result(1, 1) * result(1, 1) +
374 result(1, 2) * result(1, 2));
375 result(1, 0) /= norm1;
376 result(1, 1) /= norm1;
377 result(1, 2) /= norm1;
379 return sqrt(dxyzdXbar(0) * dxyzdXbar(0) + dxyzdXbar(1) * dxyzdXbar(1) +
380 dxyzdXbar(2) * dxyzdXbar(2));
395 dxyzdX(0, 0) += gSX(j) * nodesXYZ(j, 0);
396 dxyzdX(0, 1) += gSX(j) * nodesXYZ(j, 1);
397 dxyzdX(0, 2) += gSX(j) * nodesXYZ(j, 2);
398 dxyzdY(0, 0) += gSY(j) * nodesXYZ(j, 0);
399 dxyzdY(0, 1) += gSY(j) * nodesXYZ(j, 1);
400 dxyzdY(0, 2) += gSY(j) * nodesXYZ(j, 2);
403 result(0, 2) = dxyzdX(0, 0) * dxyzdY(0, 1) - dxyzdX(0, 1) * dxyzdY(0, 0);
404 result(0, 1) = -dxyzdX(0, 0) * dxyzdY(0, 2) + dxyzdX(0, 2) * dxyzdY(0, 0);
405 result(0, 0) = dxyzdX(0, 1) * dxyzdY(0, 2) - dxyzdX(0, 2) * dxyzdY(0, 1);
407 sqrt(result(0, 0) * result(0, 0) + result(0, 1) * result(0, 1) +
408 result(0, 2) * result(0, 2));
409 const double invNorm0 = 1. / norm0;
410 result(0, 0) *= invNorm0;
411 result(0, 1) *= invNorm0;
412 result(0, 2) *= invNorm0;
429 dxyzdX(0, 0) += gSX(j) * nodesXYZ(j, 0);
430 dxyzdX(0, 1) += gSX(j) * nodesXYZ(j, 1);
431 dxyzdX(0, 2) += gSX(j) * nodesXYZ(j, 2);
432 dxyzdY(0, 0) += gSY(j) * nodesXYZ(j, 0);
433 dxyzdY(0, 1) += gSY(j) * nodesXYZ(j, 1);
434 dxyzdY(0, 2) += gSY(j) * nodesXYZ(j, 2);
435 dxyzdZ(0, 0) += gSZ(j) * nodesXYZ(j, 0);
436 dxyzdZ(0, 1) += gSZ(j) * nodesXYZ(j, 1);
437 dxyzdZ(0, 2) += gSZ(j) * nodesXYZ(j, 2);
439 double dJ = fabs(calcDet3D(dxyzdX(0, 0), dxyzdY(0, 0), dxyzdZ(0, 0),
440 dxyzdX(0, 1), dxyzdY(0, 1), dxyzdZ(0, 1),
441 dxyzdX(0, 2), dxyzdY(0, 2), dxyzdZ(0, 2)));
456 for(
int i = 0; i < nSamplingPnts; i++) jacobian(i) = 1.;
461 if(geomNormals) normals.
setAll(*geomNormals);
464 for(
int i = 0; i < nSamplingPnts; i++) jacobian(i) = 0;
467 const double scale = 1. / invScale;
469 normals(0, 0) *=
scale;
470 normals(0, 1) *=
scale;
471 normals(0, 2) *=
scale;
474 dSMat_dX.
mult(nodesXYZ, dxyzdX);
475 for(
int i = 0; i < nSamplingPnts; i++) {
476 const double &dxdX = dxyzdX(i, 0), &dydX = dxyzdX(i, 1),
477 &dzdX = dxyzdX(i, 2);
478 const double &dxdY = normals(0, 0), &dydY = normals(0, 1),
479 &dzdY = normals(0, 2);
480 const double &dxdZ = normals(1, 0), &dydZ = normals(1, 1),
481 &dzdZ = normals(1, 2);
483 calcDet3D(dxdX, dxdY, dxdZ, dydX, dydY, dydZ, dzdX, dzdY, dzdZ);
489 if(geomNormals) normal.
setAll(*geomNormals);
492 for(
int i = 0; i < nSamplingPnts; i++) jacobian(i) = 0;
495 const double scale = 1. / invScale;
497 normal(0, 0) *=
scale;
498 normal(0, 1) *=
scale;
499 normal(0, 2) *=
scale;
502 dSMat_dX.
mult(nodesXYZ, dxyzdX);
503 dSMat_dY.
mult(nodesXYZ, dxyzdY);
504 for(
int i = 0; i < nSamplingPnts; i++) {
505 const double &dxdX = dxyzdX(i, 0), &dydX = dxyzdX(i, 1),
506 &dzdX = dxyzdX(i, 2);
507 const double &dxdY = dxyzdY(i, 0), &dydY = dxyzdY(i, 1),
508 &dzdY = dxyzdY(i, 2);
509 const double &dxdZ = normal(0, 0), &dydZ = normal(0, 1),
510 &dzdZ = normal(0, 2);
512 calcDet3D(dxdX, dxdY, dxdZ, dydX, dydY, dydZ, dzdX, dzdY, dzdZ);
518 dxyzdZ(nSamplingPnts, 3);
519 dSMat_dX.
mult(nodesXYZ, dxyzdX);
520 dSMat_dY.
mult(nodesXYZ, dxyzdY);
521 dSMat_dZ.
mult(nodesXYZ, dxyzdZ);
522 for(
int i = 0; i < nSamplingPnts; i++) {
523 const double &dxdX = dxyzdX(i, 0), &dydX = dxyzdX(i, 1),
524 &dzdX = dxyzdX(i, 2);
525 const double &dxdY = dxyzdY(i, 0), &dydY = dxyzdY(i, 1),
526 &dzdY = dxyzdY(i, 2);
527 const double &dxdZ = dxyzdZ(i, 0), &dydZ = dxyzdZ(i, 1),
528 &dzdZ = dxyzdZ(i, 2);
530 calcDet3D(dxdX, dxdY, dxdZ, dydX, dydY, dydZ, dzdX, dzdY, dzdZ);
553 const int numEl = nodesX.
size2();
554 for(
int iEl = 0; iEl < numEl; iEl++)
555 for(
int i = 0; i < nSamplingPnts; i++) jacobian(i, iEl) = 1.;
558 const int numEl = nodesX.
size2();
562 dSMat_dX.
mult(nodesX, dxdX);
563 dSMat_dX.
mult(nodesY, dydX);
564 dSMat_dX.
mult(nodesZ, dzdX);
565 for(
int iEl = 0; iEl < numEl; iEl++) {
568 nodesXYZ(i, 0) = nodesX(i, iEl);
569 nodesXYZ(i, 1) = nodesY(i, iEl);
570 nodesXYZ(i, 2) = nodesZ(i, iEl);
574 if(geomNormals) normals.
setAll(*geomNormals);
577 for(
int i = 0; i < nSamplingPnts; i++) jacobian(i, iEl) = 0;
580 const double scale = 1. / invScale;
582 normals(0, 0) *=
scale;
583 normals(0, 1) *=
scale;
584 normals(0, 2) *=
scale;
586 const double &dxdY = normals(0, 0), &dydY = normals(0, 1),
587 &dzdY = normals(0, 2);
588 const double &dxdZ = normals(1, 0), &dydZ = normals(1, 1),
589 &dzdZ = normals(1, 2);
590 for(
int i = 0; i < nSamplingPnts; i++)
591 jacobian(i, iEl) = calcDet3D(dxdX(i, iEl), dxdY, dxdZ, dydX(i, iEl),
592 dydY, dydZ, dzdX(i, iEl), dzdY, dzdZ);
596 const int numEl = nodesX.
size2();
603 dSMat_dX.
mult(nodesX, dxdX);
604 dSMat_dX.
mult(nodesY, dydX);
605 dSMat_dX.
mult(nodesZ, dzdX);
606 dSMat_dY.
mult(nodesX, dxdY);
607 dSMat_dY.
mult(nodesY, dydY);
608 dSMat_dY.
mult(nodesZ, dzdY);
609 for(
int iEl = 0; iEl < numEl; iEl++) {
612 nodesXYZ(i, 0) = nodesX(i, iEl);
613 nodesXYZ(i, 1) = nodesY(i, iEl);
614 nodesXYZ(i, 2) = nodesZ(i, iEl);
618 if(geomNormals) normal.
setAll(*geomNormals);
621 for(
int i = 0; i < nSamplingPnts; i++) jacobian(i, iEl) = 0;
624 const double scale = 1. / invScale;
626 normal(0, 0) *=
scale;
627 normal(0, 1) *=
scale;
628 normal(0, 2) *=
scale;
630 const double &dxdZ = normal(0, 0), &dydZ = normal(0, 1),
631 &dzdZ = normal(0, 2);
632 for(
int i = 0; i < nSamplingPnts; i++)
634 calcDet3D(dxdX(i, iEl), dxdY(i, iEl), dxdZ, dydX(i, iEl),
635 dydY(i, iEl), dydZ, dzdX(i, iEl), dzdY(i, iEl), dzdZ);
639 const int numEl = nodesX.
size2();
649 dSMat_dX.
mult(nodesX, dxdX);
650 dSMat_dX.
mult(nodesY, dydX);
651 dSMat_dX.
mult(nodesZ, dzdX);
652 dSMat_dY.
mult(nodesX, dxdY);
653 dSMat_dY.
mult(nodesY, dydY);
654 dSMat_dY.
mult(nodesZ, dzdY);
655 dSMat_dZ.
mult(nodesX, dxdZ);
656 dSMat_dZ.
mult(nodesY, dydZ);
657 dSMat_dZ.
mult(nodesZ, dzdZ);
658 for(
int iEl = 0; iEl < numEl; iEl++) {
659 for(
int i = 0; i < nSamplingPnts; i++)
660 jacobian(i, iEl) = calcDet3D(dxdX(i, iEl), dxdY(i, iEl), dxdZ(i, iEl),
661 dydX(i, iEl), dydY(i, iEl), dydZ(i, iEl),
662 dzdX(i, iEl), dzdY(i, iEl), dzdZ(i, iEl));
666 nodesXYZ(i, 0) = nodesX(i, iEl);
667 nodesXYZ(i, 1) = nodesY(i, iEl);
668 nodesXYZ(i, 2) = nodesZ(i, iEl);
671 for(
int i = 0; i < nSamplingPnts; i++) jacobian(i, iEl) *=
scale;
702 for(
int i = 0; i < nSamplingPnts; i++) {
713 dSMat_dX.
mult(nodesXYZ, dxyzdX);
714 for(
int i = 0; i < nSamplingPnts; i++) {
715 calcJDJ1D(dxyzdX(i, 0), normals(0, 0), normals(1, 0), dxyzdX(i, 1),
716 normals(0, 1), normals(1, 1), dxyzdX(i, 2), normals(0, 2),
723 dSMat_dX.
mult(nodesXYZ, dxyzdX);
724 dSMat_dY.
mult(nodesXYZ, dxyzdY);
725 for(
int i = 0; i < nSamplingPnts; i++) {
726 calcJDJ2D(dxyzdX(i, 0), dxyzdY(i, 0), normals(0, 0), dxyzdX(i, 1),
727 dxyzdY(i, 1), normals(0, 1), dxyzdX(i, 2), dxyzdY(i, 2),
728 normals(0, 2), i,
numMapNodes, dSMat_dX, dSMat_dY, JDJ);
735 dSMat_dX.
mult(nodesXYZ, dxyzdX);
736 dSMat_dY.
mult(nodesXYZ, dxyzdY);
737 dSMat_dZ.
mult(nodesXYZ, dxyzdZ);
738 for(
int i = 0; i < nSamplingPnts; i++) {
739 calcJDJ3D(dxyzdX(i, 0), dxyzdY(i, 0), dxyzdZ(i, 0), dxyzdX(i, 1),
740 dxyzdY(i, 1), dxyzdZ(i, 1), dxyzdX(i, 2), dxyzdY(i, 2),
741 dxyzdZ(i, 2), i,
numMapNodes, dSMat_dX, dSMat_dY, dSMat_dZ,
755 nodesXYZ, normals, JDJ);
764 SPoint3 v0(nodesXYZ(0, 0), nodesXYZ(0, 1), nodesXYZ(0, 2));
765 SPoint3 v1(nodesXYZ(1, 0), nodesXYZ(1, 1), nodesXYZ(1, 2));
766 SPoint3 v2(nodesXYZ(2, 0), nodesXYZ(2, 1), nodesXYZ(2, 2));
767 SPoint3 *IXYZ[3] = {&v0, &v1, &v2};
768 double jaci[2][2] = {
769 {IXYZ[1]->
x() - IXYZ[0]->
x(), IXYZ[2]->
x() - IXYZ[0]->
x()},
770 {IXYZ[1]->
y() - IXYZ[0]->
y(), IXYZ[2]->
y() - IXYZ[0]->
y()}};
771 double invJaci[2][2];
775 double jac[2][2] = {{0., 0.}, {0., 0.}};
779 const double dpsidx = dPhidX * invJaci[0][0] + dPhidY * invJaci[1][0];
780 const double dpsidy = dPhidX * invJaci[0][1] + dPhidY * invJaci[1][1];
781 jac[0][0] += nodesXYZ(i, 0) * dpsidx;
782 jac[0][1] += nodesXYZ(i, 0) * dpsidy;
783 jac[1][0] += nodesXYZ(i, 1) * dpsidx;
784 jac[1][1] += nodesXYZ(i, 1) * dpsidy;
786 const double dxdx = jac[0][0] * jac[0][0] + jac[0][1] * jac[0][1];
787 const double dydy = jac[1][0] * jac[1][0] + jac[1][1] * jac[1][1];
788 const double dxdy = jac[0][0] * jac[1][0] + jac[0][1] * jac[1][1];
789 const double sqr = sqrt((dxdx - dydy) * (dxdx - dydy) + 4 * dxdy * dxdy);
790 const double osqr =
sqr > 1e-8 ? 1 /
sqr : 0;
791 lambdaJ(l) = 0.5 * (dxdx + dydy -
sqr);
793 (1 - (dxdx - dydy) * osqr) * jac[0][0] - 2 * dxdy * osqr * jac[1][0];
795 (1 - (dxdx - dydy) * osqr) * jac[0][1] - 2 * dxdy * osqr * jac[1][1];
797 -2 * dxdy * osqr * jac[0][0] + (1 - (dydy - dxdx) * osqr) * jac[1][0];
799 -2 * dxdy * osqr * jac[0][1] + (1 - (dydy - dxdx) * osqr) * jac[1][1];
800 const double axixi = axx * invJaci[0][0] + axy * invJaci[0][1];
801 const double aetaeta = ayx * invJaci[1][0] + ayy * invJaci[1][1];
802 const double aetaxi = ayx * invJaci[0][0] + ayy * invJaci[0][1];
803 const double axieta = axx * invJaci[1][0] + axy * invJaci[1][1];
807 gradLambdaJ(l, i + 0 *
numMapNodes) = axixi * dPhidX + axieta * dPhidY;
808 gradLambdaJ(l, i + 1 *
numMapNodes) = aetaxi * dPhidX + aetaeta * dPhidY;
825 case TYPE_TRI:
return 2 * order - 2;
826 case TYPE_QUA:
return 2 * order - 1;
827 case TYPE_TET:
return 3 * order - 3;
828 case TYPE_PRI:
return 3 * order - 1;
829 case TYPE_HEX:
return 3 * order - 1;
831 return 3 * order - 3;
835 Msg::Error(
"Unknown element type %d, return order 0", parentType);
843 Msg::Error(
"jacobianMatrixSpace not yet implemented for pyramids");
851 case TYPE_TET: jacOrder = order - 1;
break;
854 case TYPE_HEX: jacOrder = order;
break;
856 Msg::Error(
"Unknown element type %d, return default space", type);