00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 static char const Id[] = "$Id: axesrot.c 2681 2009-07-31 04:30:53Z scholz $\n\n";
00025 static char const Td[] = "$Today: " __FILE__ " " __DATE__ " " __TIME__ " $\n\n";
00026
00027 #include "util.h"
00028 #include "petscblaslapack.h"
00029
00030 int checkortho(PetscReal *axes1)
00031 {
00032 PetscReal vec[ND];
00033
00034
00035 if (
00036 my_dnrm2(ND,axes1,1)<D_EPS ||
00037 my_dnrm2(ND,axes1+ND,1)<D_EPS ||
00038 my_dnrm2(ND,axes1+2*ND,1)<D_EPS
00039 ) {
00040 PetscPrintf(PETSC_COMM_WORLD,"rotation matrix is singular!\n");
00041 return(0);
00042 }
00043
00044
00045 if (
00046 PetscAbsReal(my_ddot(ND,axes1,1,axes1+ND,1))>D_EPS ||
00047 PetscAbsReal(my_ddot(ND,axes1,1,axes1+2*ND,1))>D_EPS ||
00048 PetscAbsReal(my_ddot(ND,axes1+ND,1,axes1+2*ND,1))>D_EPS
00049 ) {
00050 PetscPrintf(PETSC_COMM_WORLD,"rotation matrix does not form an orthogonal system!\n");
00051 return(0);
00052 }
00053
00054
00055 douter(ND,axes1,axes1+ND,vec);
00056 if (PetscAbsReal(my_ddot(ND,axes1+2*ND,1,vec,1))<D_EPS) {
00057 PetscPrintf(PETSC_COMM_WORLD,"rotation matrix left-handed or singular!\n");
00058 return(0);
00059 }
00060
00061 return(1);
00062 }
00063
00064 int AxesRotation(PetscReal *slice_p,PetscReal *slice_n,PetscReal *axes1,int inv,
00065 int n_vert,PetscReal *vertxyz,PetscReal *vertxyz2)
00066 {
00067 MagparFunctionLogBegin;
00068
00069
00070 if (slice_n!=PETSC_NULL && my_dnrm2(ND,slice_n,1)>D_EPS) {
00071
00072 PetscPrintf(PETSC_COMM_WORLD,
00073 "perpendicular plane: %g*x+%g*y+%g*z=%g\n",
00074 slice_n[0],
00075 slice_n[1],
00076 slice_n[2],
00077 my_ddot(ND,slice_n,1,slice_p,1)
00078 );
00079
00080 PetscReal t_norm;
00081 t_norm=1.0/my_dnrm2(ND,slice_n,1);
00082 PetscPrintf(PETSC_COMM_WORLD,
00083 "normalized normal: (%g,%g,%g)\n",
00084 slice_n[0]/t_norm,slice_n[1]/t_norm,slice_n[2]/t_norm
00085 );
00086
00087
00088 PetscReal axes3[ND*ND]={1,0,0,0,1,0,0,0,1};
00089 douter(ND,axes3+ND*1,slice_n,axes1);
00090 t_norm=my_dnrm2(ND,axes1,1);
00091
00092
00093 if (t_norm > 1000*D_EPS) {
00094 my_dscal(ND,1.0/t_norm,axes1,1);
00095 }
00096 else {
00097
00098 douter(ND,axes3,slice_n,axes1);
00099
00100
00101 t_norm=my_dnrm2(ND,axes1,1);
00102 assert(t_norm > D_EPS);
00103 my_dscal(ND,1.0/t_norm,axes1,1);
00104 }
00105
00106 douter(ND,slice_n,axes1,axes1+ND*1);
00107 t_norm=1.0/my_dnrm2(ND,axes1+ND*1,1);
00108 my_dscal(ND,t_norm,axes1+ND*1,1);
00109
00110 my_dcopy(ND,slice_n,1,axes1+ND*2,1);
00111 t_norm=1.0/my_dnrm2(ND,axes1+ND*2,1);
00112 my_dscal(ND,t_norm,axes1+ND*2,1);
00113
00114
00115
00116
00117 }
00118
00119
00120 if (!checkortho(axes1)) {
00121
00122
00123 douter(ND,axes1,axes1+ND,axes1+2*ND);
00124
00125 assert(checkortho(axes1));
00126 }
00127
00128 PetscPrintf(PETSC_COMM_WORLD,
00129 "rotation matrix:\n"
00130 "%10g %10g %10g\n"
00131 "%10g %10g %10g\n"
00132 "%10g %10g %10g\n",
00133 axes1[0],axes1[3],axes1[6],
00134 axes1[1],axes1[4],axes1[7],
00135 axes1[2],axes1[5],axes1[8]
00136 );
00137
00138
00139 assert(checkortho(axes1));
00140
00141
00142 ierr = PetscMemcpy(vertxyz2,vertxyz,ND*n_vert*sizeof(PetscReal));CHKERRQ(ierr);
00143
00144 if (inv==0) {
00145 PetscPrintf(PETSC_COMM_WORLD,"applying normal rotation\n");
00146
00147
00148 for (int i=0;i<n_vert;i++) {
00149 my_daxpy(ND,-1.0,slice_p,1,vertxyz2+ND*i,1);
00150 }
00151 PetscPrintf(PETSC_COMM_WORLD,
00152 "shifted by -(%g,%g,%g)\n",
00153 slice_p[0],slice_p[1],slice_p[2]
00154 );
00155
00156
00157 for (int i=0;i<n_vert;i++) {
00158 PetscReal DZero=0.0;
00159 PetscReal DOne=1.0;
00160 int IOne=1;
00161 int nd=ND;
00162 PetscReal p[ND];
00163
00164
00165
00166
00167 BLASgemv_("T",&nd,&nd,&DOne,axes1,&nd,vertxyz2+ND*i,&IOne,&DZero,p,&IOne);
00168 my_dcopy(ND,p,1,vertxyz2+ND*i,1);
00169 }
00170 }
00171 else {
00172 PetscPrintf(PETSC_COMM_WORLD,"applying inverse rotation\n");
00173
00174
00175 for (int i=0;i<n_vert;i++) {
00176 PetscReal DZero=0.0;
00177 PetscReal DOne=1.0;
00178 int IOne=1;
00179 int nd=ND;
00180 PetscReal p[ND];
00181
00182
00183
00184
00185 BLASgemv_("N",&nd,&nd,&DOne,axes1,&nd,vertxyz2+ND*i,&IOne,&DZero,p,&IOne);
00186 my_dcopy(ND,p,1,vertxyz2+ND*i,1);
00187 }
00188
00189
00190 for (int i=0;i<n_vert;i++) {
00191 my_daxpy(ND,1.0,slice_p,1,vertxyz2+ND*i,1);
00192 }
00193 PetscPrintf(PETSC_COMM_WORLD,
00194 "shifted by +(%g,%g,%g)\n",
00195 slice_p[0],slice_p[1],slice_p[2]
00196 );
00197 }
00198
00199 MagparFunctionLogReturn(0);
00200 }
00201