v0.15.0
Loading...
Searching...
No Matches
Functions
SurfacePressureComplexForLazy.cpp File Reference
#include <MoFEM.hpp>
#include <MethodForForceScaling.hpp>
#include <SurfacePressure.hpp>
#include <SurfacePressureComplexForLazy.hpp>

Go to the source code of this file.

Functions

MoFEMErrorCode Traction_hierarchical (int order, int *order_edge, double *N, double *N_face, double *N_edge[], double *t, double *t_edge[], double *t_face, double *traction, int gg)
 
MoFEMErrorCode Fext_h_hierarchical (int order, int *order_edge, double *N, double *N_face, double *N_edge[], double *diffN, double *diffN_face, double *diffN_edge[], double *t, double *t_edge[], double *t_face, double *dofs_x, double *dofs_x_edge[], double *dofs_x_face, double *idofs_x, double *idofs_x_edge[], double *idofs_x_face, double *Fext, double *Fext_egde[], double *Fext_face, double *iFext, double *iFext_egde[], double *iFext_face, int g_dim, const double *g_w)
 
MoFEMErrorCode KExt_hh_hierarchical (double eps, int order, int *order_edge, double *N, double *N_face, double *N_edge[], double *diffN, double *diffN_face, double *diffN_edge[], double *t, double *t_edge[], double *t_face, double *dofs_x, double *dofs_x_edge[], double *dofs_x_face, double *KExt_hh, double *KExt_egdeh[3], double *KExt_faceh, int g_dim, const double *g_w)
 
MoFEMErrorCode KExt_hh_hierarchical_edge (double eps, int order, int *order_edge, double *N, double *N_face, double *N_edge[], double *diffN, double *diffN_face, double *diffN_edge[], double *t, double *t_edge[], double *t_face, double *dofs_x, double *dofs_x_edge[], double *dofs_x_face, double *Khext_edge[3], double *KExt_edgeegde[3][3], double *KExt_faceedge[3], int g_dim, const double *g_w)
 
MoFEMErrorCode KExt_hh_hierarchical_face (double eps, int order, int *order_edge, double *N, double *N_face, double *N_edge[], double *diffN, double *diffN_face, double *diffN_edge[], double *t, double *t_edge[], double *t_face, double *dofs_x, double *dofs_x_edge[], double *dofs_x_face, double *KExt_hface, double *KExt_egdeface[3], double *KExt_faceface, int g_dim, const double *g_w)
 
void tetcircumcenter_tp (double a[3], double b[3], double c[3], double d[3], double circumcenter[3], double *xi, double *eta, double *zeta)
 
void tricircumcenter3d_tp (double a[3], double b[3], double c[3], double circumcenter[3], double *xi, double *eta)
 
MoFEMErrorCode KExt_hh_hierarchical (double eps, int order, int *order_edge, double *N, double *N_face, double *N_edge[], double *diffN, double *diffN_face, double *diffN_edge[], double *t, double *t_edge[], double *t_face, double *dofs_x, double *dofs_x_edge[], double *dofs_x_face, double *KExt_hh, double *KExt_edgeh[], double *KExt_faceh, int g_dim, const double *g_w)
 

Function Documentation

◆ Fext_h_hierarchical()

MoFEMErrorCode Fext_h_hierarchical ( int order,
int * order_edge,
double * N,
double * N_face,
double * N_edge[],
double * diffN,
double * diffN_face,
double * diffN_edge[],
double * t,
double * t_edge[],
double * t_face,
double * dofs_x,
double * dofs_x_edge[],
double * dofs_x_face,
double * idofs_x,
double * idofs_x_edge[],
double * idofs_x_face,
double * Fext,
double * Fext_egde[],
double * Fext_face,
double * iFext,
double * iFext_egde[],
double * iFext_face,
int g_dim,
const double * g_w )

Definition at line 618 of file SurfacePressureComplexForLazy.cpp.

625 {
626 MoFEMFunctionBeginHot;
627 int dd, nn, ee, gg;
628 if (Fext != NULL)
629 bzero(Fext, 9 * sizeof(double));
630 if (iFext != NULL)
631 bzero(iFext, 9 * sizeof(double));
632 ee = 0;
633 for (; ee < 3; ee++) {
634 int nb_dofs_edge = NBEDGE_H1(order_edge[ee]);
635 if (nb_dofs_edge == 0)
636 continue;
637 if (Fext_edge != NULL)
638 bzero(Fext_edge[ee], 3 * nb_dofs_edge * sizeof(double));
639 if (iFext_edge != NULL)
640 bzero(iFext_edge[ee], 3 * nb_dofs_edge * sizeof(double));
641 }
642 int nb_dofs_face = NBFACETRI_H1(order);
643 if (nb_dofs_face != 0) {
644 if (Fext_face != NULL)
645 bzero(Fext_face, 3 * nb_dofs_face * sizeof(double));
646 if (iFext_face != NULL)
647 bzero(iFext_face, 3 * nb_dofs_face * sizeof(double));
648 }
649 gg = 0;
650 for (; gg < g_dim; gg++) {
651 double traction[3] = {0, 0, 0};
652 CHKERR Traction_hierarchical(order, order_edge, N, N_face, N_edge, t,
653 t_edge, t_face, traction, gg);
654 __CLPK_doublecomplex xnormal[3], xs1[3], xs2[3];
655 CHKERR Normal_hierarchical(
656 // FIXME: order of tractions approximation could be different for field
657 // approx.
658 order, order_edge, order, order_edge, diffN, diffN_face, diffN_edge,
659 dofs_x, dofs_x_edge, dofs_x_face, idofs_x, idofs_x_edge, idofs_x_face,
660 xnormal, xs1, xs2, gg);
661 CHKERR Base_scale(xnormal, xs1, xs2);
662 double normal_real[3], s1_real[3], s2_real[3];
663 double normal_imag[3], s1_imag[3], s2_imag[3];
664 for (dd = 0; dd < 3; dd++) {
665 normal_real[dd] = 0.5 * xnormal[dd].r;
666 normal_imag[dd] = 0.5 * xnormal[dd].i;
667 s1_real[dd] = 0.5 * xs1[dd].r;
668 s1_imag[dd] = 0.5 * xs1[dd].i;
669 s2_real[dd] = 0.5 * xs2[dd].r;
670 s2_imag[dd] = 0.5 * xs2[dd].i;
671 }
672 nn = 0;
673 for (; nn < 3; nn++) {
674 if (Fext != NULL)
675 for (dd = 0; dd < 3; dd++) {
676 // fprintf(stderr,"%d %f %f %f %f %f %f\n",
677 // gg,g_w[gg],N[3*gg+nn],normal_real[dd],
678 // traction[0],traction[1],traction[2]);
679 Fext[3 * nn + dd] +=
680 g_w[gg] * N[3 * gg + nn] * normal_real[dd] * traction[2];
681 Fext[3 * nn + dd] +=
682 g_w[gg] * N[3 * gg + nn] * s1_real[dd] * traction[0];
683 Fext[3 * nn + dd] +=
684 g_w[gg] * N[3 * gg + nn] * s2_real[dd] * traction[1];
685 }
686 if (iFext != NULL)
687 for (dd = 0; dd < 3; dd++) {
688 iFext[3 * nn + dd] +=
689 g_w[gg] * N[3 * gg + nn] * normal_imag[dd] * traction[2];
690 iFext[3 * nn + dd] +=
691 g_w[gg] * N[3 * gg + nn] * s1_imag[dd] * traction[0];
692 iFext[3 * nn + dd] +=
693 g_w[gg] * N[3 * gg + nn] * s2_imag[dd] * traction[1];
694 }
695 }
696 ee = 0;
697 for (; ee < 3; ee++) {
698 int nb_dofs_edge = NBEDGE_H1(order_edge[ee]);
699 if (nb_dofs_edge == 0)
700 continue;
701 int nn = 0;
702 for (; nn < nb_dofs_edge; nn++) {
703 if (Fext_edge != NULL)
704 for (dd = 0; dd < 3; dd++) {
705 Fext_edge[ee][3 * nn + dd] += g_w[gg] *
706 N_edge[ee][gg * nb_dofs_edge + nn] *
707 normal_real[dd] * traction[2];
708 Fext_edge[ee][3 * nn + dd] += g_w[gg] *
709 N_edge[ee][gg * nb_dofs_edge + nn] *
710 s1_real[dd] * traction[0];
711 Fext_edge[ee][3 * nn + dd] += g_w[gg] *
712 N_edge[ee][gg * nb_dofs_edge + nn] *
713 s2_real[dd] * traction[1];
714 }
715 if (iFext_edge != NULL) {
716 for (dd = 0; dd < 3; dd++) {
717 iFext_edge[ee][3 * nn + dd] += g_w[gg] *
718 N_edge[ee][gg * nb_dofs_edge + nn] *
719 normal_imag[dd] * traction[2];
720 iFext_edge[ee][3 * nn + dd] += g_w[gg] *
721 N_edge[ee][gg * nb_dofs_edge + nn] *
722 s1_imag[dd] * traction[0];
723 iFext_edge[ee][3 * nn + dd] += g_w[gg] *
724 N_edge[ee][gg * nb_dofs_edge + nn] *
725 s2_imag[dd] * traction[1];
726 }
727 }
728 }
729 }
730 if (nb_dofs_face != 0) {
731 nn = 0;
732 for (; nn < nb_dofs_face; nn++) {
733 if (Fext_face != NULL)
734 for (dd = 0; dd < 3; dd++) {
735 Fext_face[3 * nn + dd] += g_w[gg] * N_face[gg * nb_dofs_face + nn] *
736 normal_real[dd] * traction[2];
737 Fext_face[3 * nn + dd] += g_w[gg] * N_face[gg * nb_dofs_face + nn] *
738 s1_real[dd] * traction[0];
739 Fext_face[3 * nn + dd] += g_w[gg] * N_face[gg * nb_dofs_face + nn] *
740 s2_real[dd] * traction[1];
741 }
742 if (iFext_face != NULL)
743 for (dd = 0; dd < 3; dd++) {
744 iFext_face[3 * nn + dd] += g_w[gg] *
745 N_face[gg * nb_dofs_face + nn] *
746 normal_imag[dd] * traction[2];
747 iFext_face[3 * nn + dd] += g_w[gg] *
748 N_face[gg * nb_dofs_face + nn] *
749 s1_imag[dd] * traction[0];
750 iFext_face[3 * nn + dd] += g_w[gg] *
751 N_face[gg * nb_dofs_face + nn] *
752 s2_imag[dd] * traction[1];
753 }
754 }
755 }
756 }
757 MoFEMFunctionReturnHot(0);
758}
Traction_hierarchical
MoFEMErrorCode Traction_hierarchical(int order, int *order_edge, double *N, double *N_face, double *N_edge[], double *t, double *t_edge[], double *t_face, double *traction, int gg)
Definition SurfacePressureComplexForLazy.cpp:585
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:463
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:456
order
constexpr int order
Definition dg_projection.cpp:18
Base_scale
PetscErrorCode Base_scale(__CLPK_doublecomplex *xnormal, __CLPK_doublecomplex *xs1, __CLPK_doublecomplex *xs2)
Definition fem_tools.c:741
Normal_hierarchical
PetscErrorCode Normal_hierarchical(int order_approx, int *order_edge_approx, int order, int *order_edge, double *diffN, double *diffN_face, double *diffN_edge[], double *dofs, double *dofs_edge[], double *dofs_face, double *idofs, double *idofs_edge[], double *idofs_face, __CLPK_doublecomplex *xnormal, __CLPK_doublecomplex *s1, __CLPK_doublecomplex *s2, int gg)
Complex normal.
Definition fem_tools.c:569
NBEDGE_H1
#define NBEDGE_H1(P)
Number of base function on edge for H1 space.
Definition h1_hdiv_hcurl_l2.h:55
NBFACETRI_H1
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.
Definition h1_hdiv_hcurl_l2.h:60
FTensor::dd
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)
Definition ddTensor0.hpp:33
t
constexpr double t
plate stiffness
Definition plate.cpp:58
N
const int N
Definition speed_test.cpp:3
__CLPK_doublecomplex
Definition lapack_wrap.h:34
__CLPK_doublecomplex::i
__CLPK_doublereal i
Definition lapack_wrap.h:35
__CLPK_doublecomplex::r
__CLPK_doublereal r
Definition lapack_wrap.h:35

◆ KExt_hh_hierarchical() [1/2]

MoFEMErrorCode KExt_hh_hierarchical ( double eps,
int order,
int * order_edge,
double * N,
double * N_face,
double * N_edge[],
double * diffN,
double * diffN_face,
double * diffN_edge[],
double * t,
double * t_edge[],
double * t_face,
double * dofs_x,
double * dofs_x_edge[],
double * dofs_x_face,
double * KExt_hh,
double * KExt_edgeh[],
double * KExt_faceh,
int g_dim,
const double * g_w )

Definition at line 759 of file SurfacePressureComplexForLazy.cpp.

764 {
765 MoFEMFunctionBeginHot;
766 int gg, dd, ii, nn, ee;
767 bzero(KExt_hh, 9 * 9 * sizeof(double));
768 if (KExt_edgeh != NULL) {
769 for (ee = 0; ee < 3; ee++) {
770 int nb_dofs_edge = NBEDGE_H1(order_edge[ee]);
771 bzero(KExt_edgeh[ee], 9 * 3 * nb_dofs_edge * sizeof(double));
772 }
773 }
774 int nb_dofs_face = NBFACETRI_H1(order);
775 if (KExt_faceh != NULL) {
776 bzero(KExt_faceh, 9 * 3 * nb_dofs_face * sizeof(double));
777 }
778 for (gg = 0; gg < g_dim; gg++) {
779 double traction[3] = {0, 0, 0};
780 CHKERR Traction_hierarchical(order, order_edge, N, N_face, N_edge, t,
781 t_edge, t_face, traction, gg);
782 //
783 __CLPK_doublecomplex xnormal[3], xs1[3], xs2[3];
784 double idofs_x[9];
785 for (ii = 0; ii < 9; ii++) {
786 bzero(idofs_x, 9 * sizeof(double));
787 idofs_x[ii] = eps;
788 CHKERR Normal_hierarchical(order, order_edge, // FIXME
789 order, order_edge, diffN, diffN_face,
790 diffN_edge, dofs_x, dofs_x_edge, dofs_x_face,
791 idofs_x, NULL, NULL, xnormal, xs1, xs2, gg);
792 CHKERR Base_scale(xnormal, xs1, xs2);
793 double normal_imag[3], s1_imag[3], s2_imag[3];
794 for (dd = 0; dd < 3; dd++) {
795 normal_imag[dd] = 0.5 * xnormal[dd].i / eps;
796 s1_imag[dd] = 0.5 * xs1[dd].i / eps;
797 s2_imag[dd] = 0.5 * xs2[dd].i / eps;
798 }
799 nn = 0;
800 for (; nn < 3; nn++) {
801 for (dd = 0; dd < 3; dd++) {
802 KExt_hh[ii + 9 * 3 * nn + 9 * dd] +=
803 g_w[gg] * N[3 * gg + nn] * normal_imag[dd] * traction[2];
804 KExt_hh[ii + 9 * 3 * nn + 9 * dd] +=
805 g_w[gg] * N[3 * gg + nn] * s1_imag[dd] * traction[0];
806 KExt_hh[ii + 9 * 3 * nn + 9 * dd] +=
807 g_w[gg] * N[3 * gg + nn] * s2_imag[dd] * traction[1];
808 }
809 }
810 if (KExt_edgeh != NULL) {
811 for (ee = 0; ee < 3; ee++) {
812 int nb_dofs_edge = NBEDGE_H1(order_edge[ee]);
813 for (nn = 0; nn < nb_dofs_edge; nn++) {
814 for (dd = 0; dd < 3; dd++) {
815 KExt_edgeh[ee][ii + 9 * 3 * nn + 9 * dd] +=
816 g_w[gg] * N_edge[ee][nb_dofs_edge * gg + nn] *
817 normal_imag[dd] * traction[2];
818 KExt_edgeh[ee][ii + 9 * 3 * nn + 9 * dd] +=
819 g_w[gg] * N_edge[ee][nb_dofs_edge * gg + nn] * s1_imag[dd] *
820 traction[0];
821 KExt_edgeh[ee][ii + 9 * 3 * nn + 9 * dd] +=
822 g_w[gg] * N_edge[ee][nb_dofs_edge * gg + nn] * s2_imag[dd] *
823 traction[1];
824 }
825 }
826 }
827 }
828 if (KExt_faceh != NULL) {
829 for (nn = 0; nn < nb_dofs_face; nn++) {
830 for (dd = 0; dd < 3; dd++) {
831 KExt_faceh[ii + 3 * 9 * nn + 9 * dd] +=
832 g_w[gg] * N_face[nb_dofs_face * gg + nn] * normal_imag[dd] *
833 traction[2];
834 KExt_faceh[ii + 3 * 9 * nn + 9 * dd] +=
835 g_w[gg] * N_face[nb_dofs_face * gg + nn] * s1_imag[dd] *
836 traction[0];
837 KExt_faceh[ii + 3 * 9 * nn + 9 * dd] +=
838 g_w[gg] * N_face[nb_dofs_face * gg + nn] * s2_imag[dd] *
839 traction[1];
840 }
841 }
842 }
843 }
844 }
845 MoFEMFunctionReturnHot(0);
846}
eps
static const double eps
Definition check_base_functions_derivatives_on_tet.cpp:11

◆ KExt_hh_hierarchical() [2/2]

MoFEMErrorCode KExt_hh_hierarchical ( double eps,
int order,
int * order_edge,
double * N,
double * N_face,
double * N_edge[],
double * diffN,
double * diffN_face,
double * diffN_edge[],
double * t,
double * t_edge[],
double * t_face,
double * dofs_x,
double * dofs_x_edge[],
double * dofs_x_face,
double * KExt_hh,
double * KExt_egdeh[3],
double * KExt_faceh,
int g_dim,
const double * g_w )

◆ KExt_hh_hierarchical_edge()

MoFEMErrorCode KExt_hh_hierarchical_edge ( double eps,
int order,
int * order_edge,
double * N,
double * N_face,
double * N_edge[],
double * diffN,
double * diffN_face,
double * diffN_edge[],
double * t,
double * t_edge[],
double * t_face,
double * dofs_x,
double * dofs_x_edge[],
double * dofs_x_face,
double * Khext_edge[3],
double * KExt_edgeegde[3][3],
double * KExt_faceedge[3],
int g_dim,
const double * g_w )

Definition at line 847 of file SurfacePressureComplexForLazy.cpp.

853 {
854 MoFEMFunctionBeginHot;
855 int gg, dd, ii, nn, ee, EE;
856 int nb_dofs_face = NBFACETRI_H1(order);
857 for (EE = 0; EE < 3; EE++) {
858 int nb_dofs_edge_EE = NBEDGE_H1(order_edge[EE]);
859 bzero(KExt_hedge[EE], 9 * 3 * nb_dofs_edge_EE * sizeof(double));
860 if (KExt_edgeedge != NULL) {
861 for (ee = 0; ee < 3; ee++) {
862 int nb_dofs_edge = NBEDGE_H1(order_edge[ee]);
863 bzero(KExt_edgeedge[EE][ee],
864 3 * nb_dofs_edge_EE * 3 * nb_dofs_edge * sizeof(double));
865 }
866 }
867 if (KExt_faceedge != NULL) {
868 bzero(KExt_faceedge[EE],
869 3 * nb_dofs_edge_EE * 3 * nb_dofs_face * sizeof(double));
870 }
871 }
872 for (gg = 0; gg < g_dim; gg++) {
873 double traction[3] = {0, 0, 0};
874 CHKERR Traction_hierarchical(order, order_edge, N, N_face, N_edge, t,
875 t_edge, t_face, traction, gg);
876 for (EE = 0; EE < 3; EE++) {
877 int nb_dofs_edge_EE = NBEDGE_H1(order_edge[EE]);
878 double *idofs_x_edge[3] = {NULL, NULL, NULL};
879 double idofs_x_edge_EE[3 * nb_dofs_edge_EE];
880 idofs_x_edge[EE] = idofs_x_edge_EE;
881 for (ii = 0; ii < 3 * nb_dofs_edge_EE; ii++) {
882 bzero(idofs_x_edge_EE, 3 * nb_dofs_edge_EE * sizeof(double));
883 idofs_x_edge_EE[ii] = eps;
884 __CLPK_doublecomplex xnormal[3], xs1[3], xs2[3];
885 CHKERR Normal_hierarchical(order, order_edge, // FIXME
886 order, order_edge, diffN, diffN_face,
887 diffN_edge, dofs_x, dofs_x_edge, dofs_x_face,
888 NULL, idofs_x_edge, NULL, xnormal, xs1, xs2,
889 gg);
890 CHKERR Base_scale(xnormal, xs1, xs2);
891 double normal_imag[3], s1_imag[3], s2_imag[3];
892 for (dd = 0; dd < 3; dd++) {
893 normal_imag[dd] = 0.5 * xnormal[dd].i / eps;
894 s1_imag[dd] = 0.5 * xs1[dd].i / eps;
895 s2_imag[dd] = 0.5 * xs2[dd].i / eps;
896 }
897 for (nn = 0; nn < 3; nn++) {
898 for (dd = 0; dd < 3; dd++) {
899 KExt_hedge[EE][ii + 3 * nb_dofs_edge_EE * 3 * nn +
900 3 * nb_dofs_edge_EE * dd] +=
901 g_w[gg] * N[3 * gg + nn] * normal_imag[dd] * traction[2];
902 KExt_hedge[EE][ii + 3 * nb_dofs_edge_EE * 3 * nn +
903 3 * nb_dofs_edge_EE * dd] +=
904 g_w[gg] * N[3 * gg + nn] * s1_imag[dd] * traction[0];
905 KExt_hedge[EE][ii + 3 * nb_dofs_edge_EE * 3 * nn +
906 3 * nb_dofs_edge_EE * dd] +=
907 g_w[gg] * N[3 * gg + nn] * s2_imag[dd] * traction[1];
908 }
909 }
910 if (KExt_edgeedge != NULL) {
911 for (ee = 0; ee < 3; ee++) {
912 int nb_dofs_edge = NBEDGE_H1(order_edge[ee]);
913 for (nn = 0; nn < nb_dofs_edge; nn++) {
914 for (dd = 0; dd < 3; dd++) {
915 KExt_edgeedge[EE][ee][ii + 3 * nb_dofs_edge_EE * 3 * nn +
916 3 * nb_dofs_edge_EE * dd] +=
917 g_w[gg] * N_edge[ee][nb_dofs_edge * gg + nn] *
918 normal_imag[dd] * traction[2];
919 KExt_edgeedge[EE][ee][ii + 3 * nb_dofs_edge_EE * 3 * nn +
920 3 * nb_dofs_edge_EE * dd] +=
921 g_w[gg] * N_edge[ee][nb_dofs_edge * gg + nn] * s1_imag[dd] *
922 traction[0];
923 KExt_edgeedge[EE][ee][ii + 3 * nb_dofs_edge_EE * 3 * nn +
924 3 * nb_dofs_edge_EE * dd] +=
925 g_w[gg] * N_edge[ee][nb_dofs_edge * gg + nn] * s2_imag[dd] *
926 traction[1];
927 }
928 }
929 }
930 }
931 if (KExt_faceedge != NULL) {
932 for (nn = 0; nn < nb_dofs_face; nn++) {
933 for (dd = 0; dd < 3; dd++) {
934 KExt_faceedge[EE][ii + 3 * nb_dofs_edge_EE * 3 * nn +
935 3 * nb_dofs_edge_EE * dd] +=
936 g_w[gg] * N_face[nb_dofs_face * gg + nn] * normal_imag[dd] *
937 traction[2];
938 KExt_faceedge[EE][ii + 3 * nb_dofs_edge_EE * 3 * nn +
939 3 * nb_dofs_edge_EE * dd] +=
940 g_w[gg] * N_face[nb_dofs_face * gg + nn] * s1_imag[dd] *
941 traction[0];
942 KExt_faceedge[EE][ii + 3 * nb_dofs_edge_EE * 3 * nn +
943 3 * nb_dofs_edge_EE * dd] +=
944 g_w[gg] * N_face[nb_dofs_face * gg + nn] * s2_imag[dd] *
945 traction[1];
946 }
947 }
948 }
949 }
950 }
951 }
952 MoFEMFunctionReturnHot(0);
953}

◆ KExt_hh_hierarchical_face()

MoFEMErrorCode KExt_hh_hierarchical_face ( double eps,
int order,
int * order_edge,
double * N,
double * N_face,
double * N_edge[],
double * diffN,
double * diffN_face,
double * diffN_edge[],
double * t,
double * t_edge[],
double * t_face,
double * dofs_x,
double * dofs_x_edge[],
double * dofs_x_face,
double * KExt_hface,
double * KExt_egdeface[3],
double * KExt_faceface,
int g_dim,
const double * g_w )

Definition at line 955 of file SurfacePressureComplexForLazy.cpp.

961 {
962 MoFEMFunctionBeginHot;
963 int gg, dd, ii, nn, ee;
964 int nb_dofs_face = NBFACETRI_H1(order);
965 bzero(KExt_hface, 9 * 3 * nb_dofs_face * sizeof(double));
966 if (KExt_edgeface != NULL) {
967 for (ee = 0; ee < 3; ee++) {
968 int nb_dofs_edge = NBEDGE_H1(order_edge[ee]);
969 bzero(KExt_edgeface[ee],
970 3 * nb_dofs_face * 3 * nb_dofs_edge * sizeof(double));
971 }
972 }
973 if (KExt_faceface != NULL) {
974 bzero(KExt_faceface, 3 * nb_dofs_face * 3 * nb_dofs_face * sizeof(double));
975 }
976 for (gg = 0; gg < g_dim; gg++) {
977 double traction[3] = {0, 0, 0};
978 CHKERR Traction_hierarchical(order, order_edge, N, N_face, N_edge, t,
979 t_edge, t_face, traction, gg);
980 double idofs_x_face[3 * nb_dofs_face];
981 for (ii = 0; ii < 3 * nb_dofs_face; ii++) {
982 bzero(idofs_x_face, 3 * nb_dofs_face * sizeof(double));
983 idofs_x_face[ii] = eps;
984 __CLPK_doublecomplex xnormal[3], xs1[3], xs2[3];
985 CHKERR Normal_hierarchical(
986 order, order_edge, // FIXME
987 order, order_edge, diffN, diffN_face, diffN_edge, dofs_x, dofs_x_edge,
988 dofs_x_face, NULL, NULL, idofs_x_face, xnormal, xs1, xs2, gg);
989 CHKERR Base_scale(xnormal, xs1, xs2);
990 double normal_imag[3], s1_imag[3], s2_imag[3];
991 for (dd = 0; dd < 3; dd++) {
992 normal_imag[dd] = 0.5 * xnormal[dd].i / eps;
993 s1_imag[dd] = 0.5 * xs1[dd].i / eps;
994 s2_imag[dd] = 0.5 * xs2[dd].i / eps;
995 }
996 for (nn = 0; nn < 3; nn++) {
997 for (dd = 0; dd < 3; dd++) {
998 KExt_hface[ii + 3 * nb_dofs_face * 3 * nn + 3 * nb_dofs_face * dd] +=
999 g_w[gg] * N[3 * gg + nn] * normal_imag[dd] * traction[2];
1000 KExt_hface[ii + 3 * nb_dofs_face * 3 * nn + 3 * nb_dofs_face * dd] +=
1001 g_w[gg] * N[3 * gg + nn] * s1_imag[dd] * traction[0];
1002 KExt_hface[ii + 3 * nb_dofs_face * 3 * nn + 3 * nb_dofs_face * dd] +=
1003 g_w[gg] * N[3 * gg + nn] * s2_imag[dd] * traction[1];
1004 }
1005 }
1006 if (KExt_edgeface != NULL) {
1007 for (ee = 0; ee < 3; ee++) {
1008 int nb_dofs_edge = NBEDGE_H1(order_edge[ee]);
1009 for (nn = 0; nn < nb_dofs_edge; nn++) {
1010 for (dd = 0; dd < 3; dd++) {
1011 KExt_edgeface[ee][ii + 3 * nb_dofs_face * 3 * nn +
1012 3 * nb_dofs_face * dd] +=
1013 g_w[gg] * N_edge[ee][nb_dofs_edge * gg + nn] *
1014 normal_imag[dd] * traction[2];
1015 KExt_edgeface[ee][ii + 3 * nb_dofs_face * 3 * nn +
1016 3 * nb_dofs_face * dd] +=
1017 g_w[gg] * N_edge[ee][nb_dofs_edge * gg + nn] * s1_imag[dd] *
1018 traction[0];
1019 KExt_edgeface[ee][ii + 3 * nb_dofs_face * 3 * nn +
1020 3 * nb_dofs_face * dd] +=
1021 g_w[gg] * N_edge[ee][nb_dofs_edge * gg + nn] * s2_imag[dd] *
1022 traction[1];
1023 }
1024 }
1025 }
1026 }
1027 if (KExt_faceface != NULL) {
1028 for (nn = 0; nn < nb_dofs_face; nn++) {
1029 for (dd = 0; dd < 3; dd++) {
1030 KExt_faceface[ii + 3 * nb_dofs_face * 3 * nn +
1031 3 * nb_dofs_face * dd] +=
1032 g_w[gg] * N_face[nb_dofs_face * gg + nn] * normal_imag[dd] *
1033 traction[2];
1034 KExt_faceface[ii + 3 * nb_dofs_face * 3 * nn +
1035 3 * nb_dofs_face * dd] +=
1036 g_w[gg] * N_face[nb_dofs_face * gg + nn] * s1_imag[dd] *
1037 traction[0];
1038 KExt_faceface[ii + 3 * nb_dofs_face * 3 * nn +
1039 3 * nb_dofs_face * dd] +=
1040 g_w[gg] * N_face[nb_dofs_face * gg + nn] * s2_imag[dd] *
1041 traction[1];
1042 }
1043 }
1044 }
1045 }
1046 }
1047 MoFEMFunctionReturnHot(0);
1048}

◆ tetcircumcenter_tp()

void tetcircumcenter_tp ( double a[3],
double b[3],
double c[3],
double d[3],
double circumcenter[3],
double * xi,
double * eta,
double * zeta )

◆ Traction_hierarchical()

MoFEMErrorCode Traction_hierarchical ( int order,
int * order_edge,
double * N,
double * N_face,
double * N_edge[],
double * t,
double * t_edge[],
double * t_face,
double * traction,
int gg )

Definition at line 585 of file SurfacePressureComplexForLazy.cpp.

588 {
589 MoFEMFunctionBeginHot;
590 int dd, ee;
591 for (dd = 0; dd < 3; dd++)
592 traction[dd] = cblas_ddot(3, &N[gg * 3], 1, &t[dd], 3);
593 if (t_face != NULL) {
594 int nb_dofs_face = NBFACETRI_H1(order);
595 if (nb_dofs_face > 0) {
596 for (dd = 0; dd < 3; dd++)
597 traction[dd] += cblas_ddot(nb_dofs_face, &N_face[gg * nb_dofs_face], 1,
598 &t_face[dd], 3);
599 }
600 }
601 if (t_edge != NULL) {
602 ee = 0;
603 for (; ee < 3; ee++) {
604 if (t_edge[ee] == NULL)
605 continue;
606 int nb_dofs_edge = NBEDGE_H1(order_edge[ee]);
607 if (nb_dofs_edge > 0) {
608 for (dd = 0; dd < 3; dd++) {
609 traction[dd] +=
610 cblas_ddot(nb_dofs_edge, &(N_edge[ee][gg * nb_dofs_edge]), 1,
611 &(t_edge[ee][dd]), 3);
612 }
613 }
614 }
615 }
616 MoFEMFunctionReturnHot(0);
617}

◆ tricircumcenter3d_tp()

void tricircumcenter3d_tp ( double a[3],
double b[3],
double c[3],
double circumcenter[3],
double * xi,
double * eta )
Generated by Doxygen 1.12.0 and hosted at University of Glasgow