00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 #include "field.h"
00025 #include "util/util.h"
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038 int Bele(PetscReal* bvert,PetscReal* facv1,PetscReal* facv2,PetscReal* facv3,PetscReal* matele)
00039 {
00040 PetscReal *rr,zeta[ND],zetal;
00041 PetscReal rho1[ND],rho2[ND],rho3[ND];
00042 PetscReal rho1l,rho2l,rho3l;
00043 PetscReal s1,s2,s3;
00044 PetscReal eta1[ND],eta2[ND],eta3[ND];
00045 PetscReal eta1l,eta2l,eta3l;
00046 PetscReal xi1[ND],xi2[ND],xi3[ND];
00047 PetscReal gamma1[ND],gamma2[ND],gamma3[ND];
00048 PetscReal p[ND],a,omega;
00049
00050 matele[0]=matele[1]=matele[2]=0.0;
00051
00052
00053 my_dcopy(ND,facv1,1,rho1,1);
00054 my_dcopy(ND,facv2,1,rho2,1);
00055 my_dcopy(ND,facv3,1,rho3,1);
00056
00057
00058 my_dcopy(ND,rho2,1,xi1,1);
00059 my_daxpy(ND,-1.0,rho1,1,xi1,1);
00060 my_dcopy(ND,rho3,1,xi2,1);
00061 my_daxpy(ND,-1.0,rho2,1,xi2,1);
00062 my_dcopy(ND,rho1,1,xi3,1);
00063 my_daxpy(ND,-1.0,rho3,1,xi3,1);
00064
00065
00066 douter(ND,xi1,xi2,zeta);
00067
00068
00069 zetal=my_dnrm2(ND,zeta,1);
00070 a=0.5*zetal;
00071
00072
00073 my_dscal(ND,1.0/zetal,zeta,1);
00074
00075
00076 s1=my_dnrm2(ND,xi1,1);
00077 my_dscal(ND,1.0/s1,xi1,1);
00078 s2=my_dnrm2(ND,xi2,1);
00079 my_dscal(ND,1.0/s2,xi2,1);
00080 s3=my_dnrm2(ND,xi3,1);
00081 my_dscal(ND,1.0/s3,xi3,1);
00082
00083 douter(ND,zeta,xi1,eta1);
00084 douter(ND,zeta,xi2,eta2);
00085 douter(ND,zeta,xi3,eta3);
00086
00087 gamma1[0]=gamma3[1]=my_ddot(ND,xi2,1,xi1,1);
00088 gamma1[1]=my_ddot(ND,xi2,1,xi2,1);
00089 gamma1[2]=gamma2[1]=my_ddot(ND,xi2,1,xi3,1);
00090
00091 gamma2[0]=gamma3[2]=my_ddot(ND,xi3,1,xi1,1);
00092 gamma2[2]=my_ddot(ND,xi3,1,xi3,1);
00093
00094 gamma3[0]=my_ddot(ND,xi1,1,xi1,1);
00095
00096
00097 rr=bvert;
00098
00099 PetscReal d;
00100 PointFromPlane(rr,rho1,rho2,rho3,&d);
00101 if (PetscAbsReal(d)<D_EPS) return(0);
00102
00103
00104 my_daxpy(ND,-1.0,rr,1,rho1,1);
00105 my_daxpy(ND,-1.0,rr,1,rho2,1);
00106 my_daxpy(ND,-1.0,rr,1,rho3,1);
00107
00108
00109 zetal=my_ddot(ND,zeta,1,rho1,1);
00110
00111
00112
00113
00114 if (PetscAbsReal(zetal)<=D_EPS) {
00115 return(0);
00116 }
00117
00118 rho1l=my_dnrm2(ND,rho1,1);
00119 rho2l=my_dnrm2(ND,rho2,1);
00120 rho3l=my_dnrm2(ND,rho3,1);
00121
00122 PetscReal t_nom,t_denom;
00123 t_nom=
00124 rho1l*rho2l*rho3l+
00125 rho1l*my_ddot(ND,rho2,1,rho3,1)+
00126 rho2l*my_ddot(ND,rho3,1,rho1,1)+
00127 rho3l*my_ddot(ND,rho1,1,rho2,1);
00128 t_denom=
00129 sqrt(2.0*
00130 (rho2l*rho3l+my_ddot(ND,rho2,1,rho3,1))*
00131 (rho3l*rho1l+my_ddot(ND,rho3,1,rho1,1))*
00132 (rho1l*rho2l+my_ddot(ND,rho1,1,rho2,1))
00133 );
00134
00135
00136
00137
00138
00139
00140
00141
00142 if (t_nom/t_denom<-1.0) {
00143 omega=(zetal >= 0.0 ? 1.0 : -1.0)*2.0*PETSC_PI;
00144 }
00145
00146 else if (t_nom/t_denom>1.0) {
00147 return(0);
00148 }
00149 else {
00150 omega=(zetal >= 0.0 ? 1.0 : -1.0)*2.0*acos(t_nom/t_denom);
00151 }
00152
00153 eta1l=my_ddot(ND,eta1,1,rho1,1);
00154 eta2l=my_ddot(ND,eta2,1,rho2,1);
00155 eta3l=my_ddot(ND,eta3,1,rho3,1);
00156
00157 p[0]=log((rho1l+rho2l+s1)/(rho1l+rho2l-s1));
00158 p[1]=log((rho2l+rho3l+s2)/(rho2l+rho3l-s2));
00159 p[2]=log((rho3l+rho1l+s3)/(rho3l+rho1l-s3));
00160
00161 matele[0]=(eta2l*omega-zetal*my_ddot(ND,gamma1,1,p,1))*s2/(8.0*PETSC_PI*a);
00162 matele[1]=(eta3l*omega-zetal*my_ddot(ND,gamma2,1,p,1))*s3/(8.0*PETSC_PI*a);
00163 matele[2]=(eta1l*omega-zetal*my_ddot(ND,gamma3,1,p,1))*s1/(8.0*PETSC_PI*a);
00164
00165 return(0);
00166 }