32 static PetscErrorCode
ierr;
43 det_jac *= jac[3 *
i +
i];
47 if ((
j - (
j / 2) * 2) != 0)
59 SETERRQ1(PETSC_COMM_SELF, 1,
"info = %d", info);
62 SETERRQ1(PETSC_COMM_SELF, 1,
"info = %d", info);
90 ierr = PetscMalloc(point_num *
sizeof(
double), &
w);
92 ierr = PetscMalloc(dim_num * point_num *
sizeof(
double), &x);
97 for (point = 0; point < point_num; point++) {
98 G_TRI_X[point] = x[0 + point * dim_num];
99 G_TRI_W[point] =
w[point];
134 ierr = PetscMalloc(point_num *
sizeof(
double), &
w);
136 ierr = PetscMalloc(dim_num * point_num *
sizeof(
double), &x);
141 for (point = 0; point < point_num; point++) {
142 G_TRI_X[point] = x[0 + point * dim_num];
143 G_TRI_Y[point] = x[1 + point * dim_num];
144 G_TRI_W[point] =
w[point];
174 ierr = PetscMalloc(point_num *
sizeof(
double), &
w);
176 ierr = PetscMalloc(dim_num * point_num *
sizeof(
double), &x);
180 for (point = 0; point < point_num; point++) {
181 G_TET_X[point] = x[0 + point * dim_num];
182 G_TET_Y[point] = x[1 + point * dim_num];
183 G_TET_Z[point] = x[2 + point * dim_num];
184 G_TET_W[point] =
w[point];
194 PetscErrorCode
ShapeMBTRI(
double *
N,
const double *X,
const double *Y,
198 for (; ii < G_DIM; ii++) {
199 double x = X[ii], y = Y[ii];
217 double *
normal,
double *s1,
double *s2) {
223 for (; ii < 3; ii++) {
224 diffX_ksi[ii] =
cblas_ddot(3, &coords[ii], 3, &diffN[0], 2);
225 diffX_eta[ii] =
cblas_ddot(3, &coords[ii], 3, &diffN[1], 2);
233 double Spin_diffX_ksi[9];
234 ierr =
Spin(Spin_diffX_ksi, diffX_ksi);
237 diffX_eta, 1, 0.,
normal, 1);
250 double *diff_normal) {
257 for (; ii < 3; ii++) {
258 diffX_ksi[ii] =
cblas_ddot(3, &coords[ii], 3, &diffN[0], 2);
259 diffX_eta[ii] =
cblas_ddot(3, &coords[ii], 3, &diffN[1], 2);
261 double Spin_diffX_ksi[9];
262 ierr =
Spin(Spin_diffX_ksi, diffX_ksi);
264 double Spin_diffX_eta[9];
265 ierr =
Spin(Spin_diffX_eta, diffX_eta);
268 bzero(B_ksi, 3 * 9 *
sizeof(
double));
270 bzero(B_eta, 3 * 9 *
sizeof(
double));
288 for (; ii < 3; ii++) {
289 cblas_dcopy(3, &diffN[0], 2, &B_ksi[ii * 9 + ii], 3);
290 cblas_dcopy(3, &diffN[1], 2, &B_eta[ii * 9 + ii], 3);
293 Spin_diffX_ksi, 3, B_eta, 9, 0., diff_normal, 9);
295 Spin_diffX_eta, 3, B_ksi, 9, 1., diff_normal, 9);
300 PetscErrorCode
ShapeJacMBTET(
double *diffN,
const double *coords,
double *jac) {
303 bzero(jac,
sizeof(
double) * 9);
304 for (ii = 0; ii < 4; ii++)
305 for (jj = 0; jj < 3; jj++)
306 for (kk = 0; kk < 3; kk++)
307 jac[jj * 3 + kk] += diffN[ii * 3 + kk] * coords[ii * 3 + jj];
318 PetscErrorCode
ShapeMBTET(
double *
N,
const double *G_X,
const double *G_Y,
319 const double *G_Z,
int DIM) {
322 for (; ii < DIM; ii++) {
323 double x = G_X[ii], y = G_Y[ii], z = G_Z[ii];
348 const double *elem_coords,
349 const double *glob_coords,
350 double *loc_coords) {
357 cblas_ddot(4, &diffN[0 * 3 + 0], 3, &elem_coords[0 * 3 + 0], 3);
359 cblas_ddot(4, &diffN[0 * 3 + 1], 3, &elem_coords[0 * 3 + 0], 3);
361 cblas_ddot(4, &diffN[0 * 3 + 2], 3, &elem_coords[0 * 3 + 0], 3);
363 glob_coords[0] -
cblas_ddot(4, &
N[0], 1, &elem_coords[0 * 3 + 0], 3);
367 cblas_ddot(4, &diffN[0 * 3 + 0], 3, &elem_coords[0 * 3 + 1], 3);
369 cblas_ddot(4, &diffN[0 * 3 + 1], 3, &elem_coords[0 * 3 + 1], 3);
371 cblas_ddot(4, &diffN[0 * 3 + 2], 3, &elem_coords[0 * 3 + 1], 3);
373 glob_coords[1] -
cblas_ddot(4, &
N[0], 1, &elem_coords[0 * 3 + 1], 3);
377 cblas_ddot(4, &diffN[0 * 3 + 0], 3, &elem_coords[0 * 3 + 2], 3);
379 cblas_ddot(4, &diffN[0 * 3 + 1], 3, &elem_coords[0 * 3 + 2], 3);
381 cblas_ddot(4, &diffN[0 * 3 + 2], 3, &elem_coords[0 * 3 + 2], 3);
383 glob_coords[2] -
cblas_ddot(4, &
N[0], 1, &elem_coords[0 * 3 + 2], 3);
388 SETERRQ1(PETSC_COMM_SELF, 1,
"info == %d", info);
393 const double *elem_coords,
394 const double *glob_coords,
395 double *loc_coords) {
400 A[0] =
cblas_ddot(3, &diffN[0], 2, &elem_coords[0], 2);
401 A[1] =
cblas_ddot(3, &diffN[1], 2, &elem_coords[0], 2);
402 loc_coords[0] = glob_coords[0] -
cblas_ddot(3, &
N[0], 1, &elem_coords[0], 2);
405 A[2] =
cblas_ddot(3, &diffN[0], 2, &elem_coords[1], 2);
406 A[3] =
cblas_ddot(3, &diffN[1], 2, &elem_coords[1], 2);
407 loc_coords[1] = glob_coords[1] -
cblas_ddot(3, &
N[0], 1, &elem_coords[1], 2);
410 double invA[2 * 2], detA;
411 detA = A[0] * A[3] - A[1] * A[2];
413 invA[0] = A[3] * detA;
414 invA[1] = -1.0 * A[1] * detA;
415 invA[2] = -1.0 * A[2] * detA;
416 invA[3] = A[0] * detA;
418 double loc_coords_new[2];
419 loc_coords_new[0] = invA[0] * loc_coords[0] + invA[1] * loc_coords[1];
420 loc_coords_new[1] = invA[2] * loc_coords[0] + invA[3] * loc_coords[1];
422 loc_coords[0] = loc_coords_new[0];
423 loc_coords[1] = loc_coords_new[1];
428 double *diffNinvJac) {
431 for (; ii < 4; ii++) {
433 1, 0., &diffNinvJac[ii * 3], 1);
440 for (; row < 3; row++)
441 for (col = 0; col < 3; col++) {
442 F[3 * row + col] =
cblas_ddot(4, &diffN[col], 3, &dofs[row], 3);
454 for (; ii < 4; ii++) {
456 for (jj = 0; jj < 3; jj++) {
457 tmp3[jj].
r = diffN[ii * 3 + jj];
460 cblas_zgemv(
CblasRowMajor, Trans, 3, 3, &tmp1,
invJac, 3, tmp3, 1, &tmp2,
461 &diffNinvJac[ii * 3], 1);
468 double complex diffX_x, diffX_y, diffX_z;
469 double complex diffY_x, diffY_y, diffY_z;
470 diffX_x = diffX_y = diffX_z = 0.;
471 diffY_x = diffY_y = diffY_z = 0.;
473 for (ii = 0; ii < 3; ii++) {
475 (xcoords[3 * ii + 0].
r +
I * xcoords[3 * ii + 0].
i) * diffN[2 * ii + 0];
477 (xcoords[3 * ii + 1].
r +
I * xcoords[3 * ii + 1].
i) * diffN[2 * ii + 0];
479 (xcoords[3 * ii + 2].
r +
I * xcoords[3 * ii + 2].
i) * diffN[2 * ii + 0];
481 (xcoords[3 * ii + 0].
r +
I * xcoords[3 * ii + 0].
i) * diffN[2 * ii + 1];
483 (xcoords[3 * ii + 1].
r +
I * xcoords[3 * ii + 1].
i) * diffN[2 * ii + 1];
485 (xcoords[3 * ii + 2].
r +
I * xcoords[3 * ii + 2].
i) * diffN[2 * ii + 1];
488 tmp = diffX_y * diffY_z - diffX_z * diffY_y;
489 xnormal[0].
r = creal(tmp);
490 xnormal[0].
i = cimag(tmp);
491 tmp = diffX_z * diffY_x - diffX_x * diffY_z;
492 xnormal[1].
r = creal(tmp);
493 xnormal[1].
i = cimag(tmp);
494 tmp = diffX_x * diffY_y - diffX_y * diffY_x;
495 xnormal[2].
r = creal(tmp);
496 xnormal[2].
i = cimag(tmp);
503 for (; ii < 3; ii++) {
504 for (jj = 0; jj < 3; jj++) {
505 xA[3 * ii + jj].
r = reA[3 * ii + jj];
506 xA[3 * ii + jj].
i = imA[3 * ii + jj];
519 SETERRQ(PETSC_COMM_SELF, 1,
"info == 0");
522 SETERRQ(PETSC_COMM_SELF, 1,
"info == 0");
530 SETERRQ(PETSC_COMM_SELF, 1,
"info == 0");
534 SETERRQ(PETSC_COMM_SELF, 1,
"info == 0");
543 SETERRQ(PETSC_COMM_SELF, 1,
"lapack_zgetrf(3,3,xF,3,IPIV) != 0");
545 double complex det = 1;
548 det *= xF[3 *
i +
i].
r +
I * xF[3 *
i +
i].
i;
549 if (IPIV[
i] !=
i + 1)
552 if ((
j - (
j / 2) * 2) != 0)
554 (*det_xF).r = creal(det);
555 (*det_xF).i = cimag(det);
558 PetscErrorCode
Spin(
double *spinOmega,
double *vecOmega) {
560 bzero(spinOmega, 9 *
sizeof(
double));
561 spinOmega[0 * 3 + 1] = -vecOmega[2];
562 spinOmega[0 * 3 + 2] = +vecOmega[1];
563 spinOmega[1 * 3 + 0] = +vecOmega[2];
564 spinOmega[1 * 3 + 2] = -vecOmega[0];
565 spinOmega[2 * 3 + 0] = -vecOmega[1];
566 spinOmega[2 * 3 + 1] = +vecOmega[0];
573 for (; ii < 3; ii++) {
574 for (jj = 0; jj < 3; jj++) {
575 xA[3 * ii + jj].
r = reA[3 * ii + jj];
576 xA[3 * ii + jj].
i = imA[3 * ii + jj];
582 int order_approx,
int *order_edge_approx,
int order,
int *order_edge,
583 double *diffN,
double *diffN_face,
double *diffN_edge[],
double *dofs,
584 double *dofs_edge[],
double *dofs_face,
double *idofs,
double *idofs_edge[],
590 double complex diffX_x_node, diffX_y_node, diffX_z_node;
591 double complex diffY_x_node, diffY_y_node, diffY_z_node;
598 if (dofs != NULL || idofs != NULL) {
600 for (; nn < 3; nn++) {
602 diffX_x_node += dofs[3 * nn + 0] * diffN[2 * nn + 0];
603 diffX_y_node += dofs[3 * nn + 1] * diffN[2 * nn + 0];
604 diffX_z_node += dofs[3 * nn + 2] * diffN[2 * nn + 0];
605 diffY_x_node += dofs[3 * nn + 0] * diffN[2 * nn + 1];
606 diffY_y_node += dofs[3 * nn + 1] * diffN[2 * nn + 1];
607 diffY_z_node += dofs[3 * nn + 2] * diffN[2 * nn + 1];
610 diffX_x_node +=
I * idofs[3 * nn + 0] * diffN[2 * nn + 0];
611 diffX_y_node +=
I * idofs[3 * nn + 1] * diffN[2 * nn + 0];
612 diffX_z_node +=
I * idofs[3 * nn + 2] * diffN[2 * nn + 0];
613 diffY_x_node +=
I * idofs[3 * nn + 0] * diffN[2 * nn + 1];
614 diffY_y_node +=
I * idofs[3 * nn + 1] * diffN[2 * nn + 1];
615 diffY_z_node +=
I * idofs[3 * nn + 2] * diffN[2 * nn + 1];
619 double complex diffX_x, diffX_y, diffX_z;
620 double complex diffY_x, diffY_y, diffY_z;
621 diffX_x = diffX_x_node;
622 diffX_y = diffX_y_node;
623 diffX_z = diffX_z_node;
624 diffY_x = diffY_x_node;
625 diffY_y = diffY_y_node;
626 diffY_z = diffY_z_node;
627 if (dofs_face != NULL || idofs_face != NULL) {
630 if (nb_dofs_face > 0) {
631 if (dofs_face != NULL) {
632 diffX_x +=
cblas_ddot(nb_dofs_face, &dofs_face[0], 3,
633 &diffN_face[gg * 2 * nb_dofs_approx_face + 0], 2);
634 diffX_y +=
cblas_ddot(nb_dofs_face, &dofs_face[1], 3,
635 &diffN_face[gg * 2 * nb_dofs_approx_face + 0], 2);
636 diffX_z +=
cblas_ddot(nb_dofs_face, &dofs_face[2], 3,
637 &diffN_face[gg * 2 * nb_dofs_approx_face + 0], 2);
638 diffY_x +=
cblas_ddot(nb_dofs_face, &dofs_face[0], 3,
639 &diffN_face[gg * 2 * nb_dofs_approx_face + 1], 2);
640 diffY_y +=
cblas_ddot(nb_dofs_face, &dofs_face[1], 3,
641 &diffN_face[gg * 2 * nb_dofs_approx_face + 1], 2);
642 diffY_z +=
cblas_ddot(nb_dofs_face, &dofs_face[2], 3,
643 &diffN_face[gg * 2 * nb_dofs_approx_face + 1], 2);
645 if (idofs_face != NULL) {
648 &diffN_face[gg * 2 * nb_dofs_approx_face + 0], 2);
651 &diffN_face[gg * 2 * nb_dofs_approx_face + 0], 2);
654 &diffN_face[gg * 2 * nb_dofs_approx_face + 0], 2);
657 &diffN_face[gg * 2 * nb_dofs_approx_face + 1], 2);
660 &diffN_face[gg * 2 * nb_dofs_approx_face + 1], 2);
663 &diffN_face[gg * 2 * nb_dofs_approx_face + 1], 2);
668 if (dofs_edge != NULL || idofs_edge != NULL) {
669 for (; ee < 3; ee++) {
670 int nb_dofs_edge =
NBEDGE_H1(order_edge[ee]);
671 int nb_dofs_approx_edge =
NBEDGE_H1(order_edge_approx[ee]);
672 if (nb_dofs_edge > 0) {
673 if (dofs_edge != NULL) {
674 if (dofs_edge[ee] != NULL) {
676 nb_dofs_edge, &(dofs_edge[ee])[0], 3,
677 &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 0], 2);
679 nb_dofs_edge, &(dofs_edge[ee])[1], 3,
680 &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 0], 2);
682 nb_dofs_edge, &(dofs_edge[ee])[2], 3,
683 &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 0], 2);
685 nb_dofs_edge, &(dofs_edge[ee])[0], 3,
686 &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 1], 2);
688 nb_dofs_edge, &(dofs_edge[ee])[1], 3,
689 &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 1], 2);
691 nb_dofs_edge, &(dofs_edge[ee])[2], 3,
692 &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 1], 2);
695 if (idofs_edge != NULL) {
696 if (idofs_edge[ee] == NULL)
700 nb_dofs_edge, &(idofs_edge[ee])[0], 3,
701 &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 0], 2);
704 nb_dofs_edge, &(idofs_edge[ee])[1], 3,
705 &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 0], 2);
708 nb_dofs_edge, &(idofs_edge[ee])[2], 3,
709 &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 0], 2);
712 nb_dofs_edge, &(idofs_edge[ee])[0], 3,
713 &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 1], 2);
716 nb_dofs_edge, &(idofs_edge[ee])[1], 3,
717 &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 1], 2);
720 nb_dofs_edge, &(idofs_edge[ee])[2], 3,
721 &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 1], 2);
727 normal[0] = diffX_y * diffY_z - diffX_z * diffY_y;
728 normal[1] = diffX_z * diffY_x - diffX_x * diffY_z;
729 normal[2] = diffX_x * diffY_y - diffX_y * diffY_x;
731 for (;
dd < 3;
dd++) {
736 xs1[0].
r = creal(diffX_x);
737 xs1[0].
i = cimag(diffX_x);
738 xs1[1].
r = creal(diffX_y);
739 xs1[1].
i = cimag(diffX_y);
740 xs1[2].
r = creal(diffX_z);
741 xs1[2].
i = cimag(diffX_z);
744 xs2[0].
r = creal(diffY_x);
745 xs2[0].
i = cimag(diffY_x);
746 xs2[1].
r = creal(diffY_y);
747 xs2[1].
i = cimag(diffY_y);
748 xs2[2].
r = creal(diffY_z);
749 xs2[2].
i = cimag(diffY_z);
757 complex
double xnrm2_normal = csqrt(cpow(xnormal[0].
r +
I * xnormal[0].
i, 2) +
758 cpow(xnormal[1].
r +
I * xnormal[1].
i, 2) +
759 cpow(xnormal[2].
r +
I * xnormal[2].
i, 2));
761 for (;
dd < 3;
dd++) {
762 complex
double s1 = (xs1[
dd].
r +
I * xs1[
dd].
i) * xnrm2_normal;
763 complex
double s2 = (xs2[
dd].
r +
I * xs2[
dd].
i) * xnrm2_normal;
764 xs1[
dd].
r = creal(s1);
765 xs1[
dd].
i = cimag(s1);
766 xs2[
dd].
r = creal(s2);
767 xs2[
dd].
i = cimag(s2);
776 for (; ii < DIM; ii++) {
792 PetscErrorCode
ShapeMBTRIQ(
double *
N,
const double *X,
const double *Y,
796 for (; ii < G_DIM; ii++) {
797 double x = X[ii], y = Y[ii];
811 for (; ii < G_DIM; ii++) {
812 double x = X[ii], y = Y[ii];
830 #define N_MBTETQ0(x, y, z) ((2. * (1. - x - y - z) - 1.) * (1. - x - y - z))
831 #define N_MBTETQ1(x, y, z) ((2. * x - 1.) * x)
832 #define N_MBTETQ2(x, y, z) ((2. * y - 1.) * y)
833 #define N_MBTETQ3(x, y, z) ((2. * z - 1.) * z)
834 #define N_MBTETQ4(x, y, z) (4. * (1. - x - y - z) * x)
835 #define N_MBTETQ5(x, y, z) (4. * x * y)
836 #define N_MBTETQ6(x, y, z) (4. * (1. - x - y - z) * y)
837 #define N_MBTETQ7(x, y, z) (4. * (1. - x - y - z) * z)
838 #define N_MBTETQ8(x, y, z) (4. * x * z)
839 #define N_MBTETQ9(x, y, z) (4. * y * z)
840 #define diffN_MBTETQ0x(x, y, z) (-3. + 4. * x + 4. * y + 4. * z)
841 #define diffN_MBTETQ0y(x, y, z) (-3. + 4. * x + 4. * y + 4. * z)
842 #define diffN_MBTETQ0z(x, y, z) (-3. + 4. * x + 4. * y + 4. * z)
843 #define diffN_MBTETQ1x(x, y, z) (4. * x - 1.)
844 #define diffN_MBTETQ1y(x, y, z) (0.)
845 #define diffN_MBTETQ1z(x, y, z) (0.)
846 #define diffN_MBTETQ2x(x, y, z) (0.)
847 #define diffN_MBTETQ2y(x, y, z) (4. * y - 1.)
848 #define diffN_MBTETQ2z(x, y, z) (0.)
849 #define diffN_MBTETQ3x(x, y, z) (0.)
850 #define diffN_MBTETQ3y(x, y, z) (0.)
851 #define diffN_MBTETQ3z(x, y, z) (4. * z - 1.)
852 #define diffN_MBTETQ4x(x, y, z) (-8. * x + 4. - 4. * y - 4. * z)
853 #define diffN_MBTETQ4y(x, y, z) (-4. * x)
854 #define diffN_MBTETQ4z(x, y, z) (-4. * x)
855 #define diffN_MBTETQ5x(x, y, z) (4. * y)
856 #define diffN_MBTETQ5y(x, y, z) (4. * x)
857 #define diffN_MBTETQ5z(x, y, z) (0.)
858 #define diffN_MBTETQ6x(x, y, z) (-4. * y)
859 #define diffN_MBTETQ6y(x, y, z) (-8. * y + 4. - 4. * x - 4. * z)
860 #define diffN_MBTETQ6z(x, y, z) (-4. * y)
861 #define diffN_MBTETQ7x(x, y, z) (-4. * z)
862 #define diffN_MBTETQ7y(x, y, z) (-4. * z)
863 #define diffN_MBTETQ7z(x, y, z) (-8. * z + 4. - 4. * x - 4. * y)
864 #define diffN_MBTETQ8x(x, y, z) (4. * z)
865 #define diffN_MBTETQ8y(x, y, z) (0.)
866 #define diffN_MBTETQ8z(x, y, z) (4. * x)
867 #define diffN_MBTETQ9x(x, y, z) (0.)
868 #define diffN_MBTETQ9y(x, y, z) (4. * z)
869 #define diffN_MBTETQ9z(x, y, z) (4. * y)
921 const double *Z,
const int G_DIM) {
924 for (; ii < G_DIM; ii++) {
925 double x = X[ii], y = Y[ii], z = Z[ii];
940 const double *Y,
const double *Z,
944 for (; ii < G_DIM; ii++) {
945 double x = X[ii], y = Y[ii], z = Z[ii];
983 bzero(Jac,
sizeof(
double) * 9);
984 for (ii = 0; ii < 10; ii++)
985 for (jj = 0; jj < 3; jj++)
986 for (kk = 0; kk < 3; kk++)
987 Jac[jj * 3 + kk] += diffN[ii * 3 + kk] * coords[ii * 3 + jj];
992 const double *diffN,
const double *coords,
998 for (; ii < G_DIM; ii++) {
1010 double detJac_at_Gauss_Points[G_DIM];
1014 for (; ii < G_DIM; ii++) {
1015 vol += G_TET_W[ii] * (detJac_at_Gauss_Points[ii]) / 6;
1020 const double *elem_coords,
1021 const double *glob_coords,
1022 double *loc_coords,
const double eps) {
1027 float NORM_dR = 1000.;
1031 R[0] = glob_coords[0] -
cblas_ddot(10, &
N[0], 1, &elem_coords[0], 3);
1032 R[1] = glob_coords[1] -
cblas_ddot(10, &
N[0], 1, &elem_coords[1], 3);
1033 R[2] = glob_coords[2] -
cblas_ddot(10, &
N[0], 1, &elem_coords[2], 3);
1035 while ((NORM_dR / NORM_R0) >
eps) {
1038 A[0 + 3 * 0] =
cblas_ddot(10, &diffN[0 * 3 + 0], 3, &elem_coords[0], 3);
1039 A[0 + 3 * 1] =
cblas_ddot(10, &diffN[0 * 3 + 1], 3, &elem_coords[0], 3);
1040 A[0 + 3 * 2] =
cblas_ddot(10, &diffN[0 * 3 + 2], 3, &elem_coords[0], 3);
1041 R[0] = glob_coords[0] -
cblas_ddot(10, &
N[0], 1, &elem_coords[0], 3);
1043 A[1 + 3 * 0] =
cblas_ddot(10, &diffN[0 * 3 + 0], 3, &elem_coords[1], 3);
1044 A[1 + 3 * 1] =
cblas_ddot(10, &diffN[0 * 3 + 1], 3, &elem_coords[1], 3);
1045 A[1 + 3 * 2] =
cblas_ddot(10, &diffN[0 * 3 + 2], 3, &elem_coords[1], 3);
1046 R[1] = glob_coords[1] -
cblas_ddot(10, &
N[0], 1, &elem_coords[1], 3);
1049 cblas_ddot(10, &diffN[0 * 3 + 0], 3, &elem_coords[0 * 3 + 2], 3);
1051 cblas_ddot(10, &diffN[0 * 3 + 1], 3, &elem_coords[0 * 3 + 2], 3);
1053 cblas_ddot(10, &diffN[0 * 3 + 2], 3, &elem_coords[0 * 3 + 2], 3);
1054 R[2] = glob_coords[2] -
cblas_ddot(10, &
N[0], 1, &elem_coords[2], 3);
1060 ShapeMBTETQ(
N, loc_coords[0], loc_coords[1], loc_coords[2]);
double cblas_ddot(const int N, const double *X, const int incX, const double *Y, const int incY)
double cblas_dnrm2(const int N, const double *X, const int incX)
void cblas_dgemv(const enum CBLAS_ORDER order, const enum CBLAS_TRANSPOSE TransA, const int M, const int N, const double alpha, const double *A, const int lda, const double *X, const int incX, const double beta, double *Y, const int incY)
void cblas_dcopy(const int N, const double *X, const int incX, double *Y, const int incY)
void cblas_daxpy(const int N, const double alpha, const double *X, const int incX, double *Y, const int incY)
void cblas_dgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, const int M, const int N, const int K, const double alpha, const double *A, const int lda, const double *B, const int ldb, const double beta, double *C, const int ldc)
void cblas_zgemv(const enum CBLAS_ORDER order, const enum CBLAS_TRANSPOSE TransA, const int M, const int N, const void *alpha, const void *A, const int lda, const void *X, const int incX, const void *beta, void *Y, const int incY)
useful compiler directives and definitions
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
void gm_rule_set(int rule, int dim_num, int point_num, double w[], double x[])
int gm_rule_size(int rule, int dim_num)
Functions to approximate hierarchical spaces.
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.
static __CLPK_integer lapack_zgetri(__CLPK_integer n, __CLPK_doublecomplex *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublecomplex *work, __CLPK_integer lwork)
static __CLPK_integer lapack_zpotri(char uplo, __CLPK_integer n, __CLPK_doublecomplex *a, __CLPK_integer lda)
static __CLPK_integer lapack_zgetrf(__CLPK_integer m, __CLPK_integer n, __CLPK_doublecomplex *a, __CLPK_integer lda, __CLPK_integer *ipiv)
static __CLPK_integer lapack_dgesv(__CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *b, __CLPK_integer ldb)
static __CLPK_integer lapack_dgetrf(__CLPK_integer m, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv)
static __CLPK_integer lapack_dgetri(__CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *work, __CLPK_integer lwork)
static __CLPK_integer lapack_zpotrf(char uplo, __CLPK_integer n, __CLPK_doublecomplex *a, __CLPK_integer lda)
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
const double r
rate factor