|
Bemsolver 2.0
|
bem.h
Go to the documentation of this file.
00001 //Todo heap stack problem with D3ImportedElectrode 00002 00003 //this is already done: 00004 //to include fortran libraries read especially "Intel Fortran/C Mixed-Language Programs Overview" in Intel Visual Fortran User Guide Vol I 00005 //but it also works without including libifcore.lib 00006 //you just have to set in the ctest Project property tab at linker/input/additional dependencies to ..\Lib2\Debug\Lib2.lib or ..\Lib2\Release\Lib2.lib 00007 //and the linker/common/additional library directory to D:\Programme\Intel\Fortran\compiler80\IA32\LIB 00008 //rightclick on Projectmappe to set dependencies and define the file to be debugged and started 00009 00010 //if linker does not find ifconsol.lib then change setting in properties of whole project/linker/allgemein/zusaetzliche bibliotheksverzeichnisse to the proper directory of the intel fortran compiler 00011 00012 // testlapack.cpp : Definiert den Einstiegspunkt für die Konsolenanwendung. 00013 //in the fortran project: 00014 //Remove "/libdir:noauto". In the project property pages, this is Fortran..Libraries..Disable default library search rules - change to "No". 00015 //This setting, which is "Yes" by default in static library projects, tells the compiler to not include in the object file directives to pull in the run-time libraries. This is correct when the library is used by another Fortran project, but not when the main program is not Fortran. 00016 //In order to avoid library collisions also Librarian..General..Ignore all default libraries to yes 00017 //Take care that the project type in Fortran..Libraries is set to multithreaded debug/ multithreaded release respectively (if the current c project has the same settings) 00018 00019 00020 //Euler angles as defined in goldstein page 146. Best understood as: 00021 //phi, theta psi: 00022 //first rotate object around z axis with angle psi 00023 //then rotate in the objects rotated frame around the local x axis with angle theta 00024 //finalay rotat in the objects rotated frame around the local z axis. 00025 //todo put spaceunits into checksum 00026 00027 #ifndef __BEM_H 00028 #define __BEM_H 00029 #include <complex> //modified std::complex with inversed storage order or real and imag in order to comply with fftw_complex 00030 using namespace std; 00031 00032 #ifdef WIN32 00033 #include <Windows4Root.h>//is needed otherwise TCanvas.h creates an error at CreateWindow, because this is a macro in MFC 00034 #endif 00035 00036 00037 #undef GetFreeSpace 00038 #include <TGeoMatrix.h> 00039 #include <TGeoManager.h> 00040 #include <TGeoArb8.h> 00041 #include "dl_creationadapter.h" 00042 #ifndef __CINT__ 00043 00044 #define BEMREVISION 2 00045 00046 00047 #include "lapack.h" 00048 #include <float.h> 00049 #include <fftw3.h> 00050 #include <gsl/gsl_sf_bessel.h> 00051 #include <gsl/gsl_sf_ellint.h> 00052 #include <math.h> 00053 #include <gsl/gsl_integration.h> 00054 #include <gsl/gsl_matrix.h> 00055 #include <gsl/gsl_vector.h> 00056 #include <gsl/gsl_linalg.h> 00057 00058 00059 #include "stdafx.h" 00060 #else 00061 typedef double mxArray; 00062 #endif 00063 00064 00065 #include <vector> 00066 #include <list> 00067 00068 //#define sqr(x) pow((x),2) 00069 #define sqr(x) ((x)*(x)) 00070 00071 const double pi=3.1415926535897932384626433832795028841971693993751058209749445923078164062862090; 00072 00073 // solving the matrix equation A*x=b using LAPACK 00074 00075 00076 #define size_ 3 // dimension of matrix 00077 00078 double testfn(double x,double y,double z); 00079 00080 void q(); 00081 00082 00083 typedef double trianglefn(double,double,double,double,double,double); 00084 00085 class calcD3triangle{ 00086 public: 00087 calcD3triangle(); 00088 ~calcD3triangle(); 00089 double GetIntL1(double x1,double y1,double z1,double x2,double y2,double z2,double x3,double y3,double z3); 00090 double GetIntL1OutOfPlane(double x1,double y1,double z1,double x2,double y2,double z2,double x3,double y3,double z3,double rp); 00091 double GetIntL1divR3(double x1,double y1,double z1,double x2,double y2,double z2,double x3,double y3,double z3,double rp); 00092 //triangleint(&testfn,0,0,0, 1,0,0, 0,1,0); 00093 00094 double triangleint(trianglefn *fn,double x0,double y0,double z0,double x1,double y1,double z1,double x2,double y2,double z2,double x3,double y3,double z3); 00095 static double G3D(double x0,double y0,double z0,double x1,double y1,double z1); 00096 static double dG3Ddx(double x0,double y0,double z0,double x1,double y1,double z1); 00097 static double dG3Ddy(double x0,double y0,double z0,double x1,double y1,double z1); 00098 static double dG3Ddz(double x0,double y0,double z0,double x1,double y1,double z1); 00099 double G3Danalytic(double xp,double yp,double zp,double x1,double y1,double z1,double x2,double y2,double z2,double x3,double y3,double z3,bool inversenorm); 00100 double G3DdnAnalytic(double xp,double yp,double zp,double x1,double y1,double z1,double x2,double y2,double z2,double x3,double y3,double z3,bool inversenorm); 00101 void triangleint_pot_exyz(double x0,double y0,double z0,double x1,double y1,double z1,double x2,double y2,double z2,double x3,double y3,double z3,bool inversenorm,double &pot,double &ex,double &ey,double &ez); 00102 00103 protected: 00104 void gauleg(double x1, double x2, double x[], double w[], int n); //from nr c p152 c4-5.pdf 00105 //Given the lower and upper limits of integration x1 and x2, and given n, this routine returns 00106 //arrays x[1..n] and w[1..n] of length n, containing the abscissas and weights of the Gauss- 00107 //Legendre n-point quadrature formula. 00108 00109 void gauslegsimplextriangle(double xi[],double eta[],double c[],int n); 00110 00111 double EPS; //EPS is the relative precision. 00112 double *xi; 00113 double *eta; 00114 double *w; 00115 int numquad; 00116 int n; 00117 }; 00118 00119 00120 00121 00122 //x father class of everything 00123 00124 //D3element needs list of sollutions for each 0 ...0 1 0...0 electrode voltage configuration 00125 //new electrode container class 00126 //new world container class 00127 //MatlabDisplay needs other depth than real element depth 00128 class D3element; 00129 typedef D3element *PD3element; 00130 typedef list<PD3element> PD3elementList; 00131 class D3element{ 00132 public: 00161 virtual void correctNorm(double x0,double y0,double z0); 00162 virtual bool IntersectWithRay(double x0,double y0,double z0,double xdir,double ydir,double zdir,double &mul); 00163 virtual void InsertTGeoVolume(TGeoVolume *top,TGeoMedium *mat,TGeoMedium *matVak,TGeoManager *geom); 00164 virtual double GetSelfPotential(); 00165 virtual double GetPotentialAt(double x,double y,double z); 00166 virtual double GetSelfDoubleLayerPotential(); 00167 virtual double GetDoubleLayerPotentialAt(double x,double y,double z); 00168 virtual void GetPotentialAndFieldAt(double x,double y,double z,double &pot,double &ex,double &ey,double &ez); 00169 virtual double GetArea(); 00170 virtual void GetCenter(double &x,double &y,double &z); 00171 D3element( double x_,double y_,double z_, 00172 double phi_,double theta_,double psi_, 00173 bool inversenorm_=false,bool refineable_=true,D3element* parent_=NULL,double epsilon=0,string name=""); 00174 ~D3element(); 00175 virtual void GetReferencePoint(double &x,double &y, double &z); 00176 virtual void init(bool ignorefirst=false); 00195 virtual void refine(double length); 00196 virtual void refine(double length,int num); //?? unterschied zu oben? -> evtl um sachen zu überspringen? //x 00197 int GetAmountOfSubelements(); 00198 int PrintAmountOfSubelements(); 00199 void GetListOfBaseElements(PD3element *el,int &cnt); 00200 virtual void rotate(double phi_,double theta_,double psi_); 00201 virtual void shift(double xs,double ys,double zs); 00202 virtual void rotate2(double phi_,double theta_,double psi_,double x,double y,double z,double &x2,double &y2,double &z2); 00203 virtual void rotateinv(double phi_,double theta_,double psi_,double x,double y,double z,double &x2,double &y2,double &z2); 00204 00205 // virtual int Get3DMatlab( mxArray *plhs[],int cnt=-1); 00206 D3element *parent; 00207 void Add(PD3element el); 00208 virtual void GetTriangle(double *A,double *B,double *C,double *COL); 00209 virtual void GetRectangle(double *A,double *B,double *C,double *D,double *COL); 00210 virtual void SetNormTowards(double x,double y,double z,bool towards); 00211 int Color; 00212 virtual void createNewSubelements(double length); 00213 virtual string getName(); 00214 virtual void setName(string name2); 00215 virtual bool GetInversenorm(); 00216 string name; 00217 list<D3element*> subelement; 00218 protected: 00219 double epsilon; 00220 TGeoHMatrix h; 00221 virtual void deleteSubelements(); 00222 double x; 00223 double y; 00224 double z; 00225 double phi; 00226 double theta; 00227 double psi; 00228 bool isBaseElement; 00229 00230 bool refineable; 00231 bool inversenorm; 00232 00233 }; 00234 00235 00236 00237 class D3electrode; 00238 00239 class D3electrode:public D3element{ 00240 public: 00246 D3electrode(); 00247 //D3electrode(const D3electrode &cp); 00248 ~D3electrode(); 00254 void insert(PD3element el); 00260 void SetVoltage(double voltage_); 00264 double GetVoltage(); 00290 void SetSymmetryWith(bool xsym,bool ysym,bool zsym,bool diagmirror,D3electrode *symel2,int rotsym=0,int totalrot=0); 00313 void SetSymmetry(bool xsym,bool ysym,bool zsym,bool diagmirror,int rotsym=0,int totalrot2=0); 00314 D3electrode *GetSymmetryWith(int num); 00315 int GetTotalrot(); 00316 00317 void SetCardinalNumber(int num); 00318 int GetCardinalNumber(); 00319 protected: 00320 int cardinalnum; 00321 double voltage; 00322 D3electrode **symel; 00323 int totalrot; 00324 }; 00325 00326 class D3sortedelement{ 00327 public: 00328 D3sortedelement():done(false){} 00329 ~D3sortedelement(){} 00330 friend bool SmallerThan(D3sortedelement &a,D3sortedelement &b,double eps){ 00331 if(abs(a.axis-b.axis)<eps) { 00332 if(a.radius<b.radius) return true; 00333 else return false; 00334 } 00335 else if(a.axis<b.axis) return true; 00336 else return false; 00337 00338 } 00339 bool done; 00340 double axis;//Major sort value 00341 double radius;//minor sort value 00342 int index; //für dfdnAll und x 00343 D3electrode *electrodeptr; //Ptr auf die das Element beinhaltende Elektrode 00344 }; 00345 00346 class D3sorter{ 00347 private: 00348 D3sortedelement *a; 00349 int n; 00350 public: 00351 void sort(D3sortedelement *arr,int size,double epsi) 00352 { 00353 a=arr; 00354 n=size; 00355 eps=epsi; 00356 quicksort(0, n-1); 00357 } 00358 00359 void quicksort (int l, int r) 00360 { 00361 if(r>l){ 00362 int i=l-1, j=r, tmp; 00363 if(r-l > 3){ //Median of three 00364 int m=l+(r-l)/2; 00365 D3sortedelement &x=a[m]; 00366 if(SmallerThan(x,a[l],eps)) exchange(m,l); 00367 if(SmallerThan(a[r],a[l],eps)) exchange(r,l); 00368 else if(SmallerThan(a[m],a[r],eps)) exchange(r,m); 00369 } 00370 00371 for(;;){ 00372 while(SmallerThan(a[++i],a[r],eps)); 00373 while(SmallerThan(a[r],a[--j],eps) && j>i); 00374 if(i>=j) break; 00375 exchange(i,j); 00376 } 00377 exchange(i,r); 00378 00379 00380 quicksort( l, i-1); 00381 quicksort( i+1, r); 00382 00383 } 00384 00385 }; 00386 00387 void exchange(int i, int j) 00388 { 00389 D3sortedelement t=a[i]; 00390 a[i]=a[j]; 00391 a[j]=t; 00392 } 00393 protected: 00394 double eps; 00395 00396 00397 }; 00398 00399 00400 00401 //x 00402 class D3ImportedElectrodes:public DL_CreationAdapter{ 00403 public: 00421 D3ImportedElectrodes(); 00422 00423 */ 00424 ~D3ImportedElectrodes(); 00435 bool Import(const char* file,bool ignore3DFace_=false,bool ignorePolyline_=false); 00436 00437 virtual void addLayer(const DL_LayerData& data); 00438 virtual void addPoint(const DL_PointData& data); 00439 virtual void addLine(const DL_LineData& data); 00440 virtual void addArc(const DL_ArcData& data); 00441 virtual void addCircle(const DL_CircleData& data); 00442 virtual void addPolyline(const DL_PolylineData& data); 00443 virtual void addVertex(const DL_VertexData& data); 00444 virtual void add3DFace(const DL_3DFaceData &data); 00445 virtual void addBlockRecord(const DL_BlockRecordData &data); 00446 void endSequence(); 00450 void ListElectrodes(); 00451 D3electrode &FindElectrode(const char* name); 00452 void printAttributes(); 00453 map<string,D3electrode> electrode; 00454 protected: 00455 int Model_Space_Handle; 00456 bool ignore_Model_Space_Handle; 00457 bool ignore3DFace; 00458 bool ignorePolyline; 00459 bool importingPolyline; 00460 int vertexcnt; 00461 double x[4]; 00462 double y[4]; 00463 double z[4]; 00464 }; 00465 // 00466 // /\ 00467 // / \ 00468 // / \ 00469 // ^ / \ 00470 // x3 /________\ 00471 // / \ / \ 00472 // / \ / \ 00473 // / \ / \ 00474 // /_______\/_______\ 00475 // xyz x2> 00476 // 00477 //phi,theta,psi are euler coordinates. rotation point is (x,y,z) 00478 00479 00480 class D3triangle:public D3element{ 00481 public: 00482 virtual void InsertTGeoVolume(TGeoVolume *top,TGeoMedium *matVak,TGeoMedium *mat,TGeoManager *geom); 00483 virtual bool IntersectWithRay(double x0,double y0,double z0,double xdir,double ydir,double zdir,double &mul); 00484 static bool IntersectWithRay(double ax,double ay,double az,double bx,double by,double bz,double cx,double cy,double cz,double xdir,double ydir,double zdir,double &mul); 00485 00486 virtual void SetNormTowards(double x,double y,double z,bool towards); 00487 virtual void GetReferencePoint(double &x,double &y, double &z); 00488 virtual double GetArea(); 00489 virtual void GetCenter(double &x,double &y,double &z); 00490 virtual double GetSelfPotential(); 00491 virtual inline double GetPotentialAt(double x,double y,double z); 00492 virtual void GetPotentialAndFieldAt(double x,double y,double z,double &pot,double &ex,double &ey,double &ez); 00493 virtual double GetSelfDoubleLayerPotential(); 00494 virtual inline double GetDoubleLayerPotentialAt(double x,double y,double z); 00513 00514 D3triangle( double xa_,double ya_,double xb_,double yb_, 00515 double x,double y,double z, 00516 double phi_,double theta_,double psi_,bool inversenorm=false,PD3element parent=NULL,bool stripeit=true); 00522 00523 D3triangle(double x1_,double y1_,double z1_,double x2_,double y2_,double z2_,double x3_,double y3_,double z3_, 00524 bool inversenorm=false,PD3element parent=NULL,bool stripeit=true); 00525 00526 ~D3triangle(); 00527 virtual void rotate(double phi_,double theta_,double psi_); 00528 virtual void shift(double xs,double ys,double zs); 00529 virtual void GetTriangle(double *A,double *B,double *C,double *COL); 00530 virtual void createNewSubelements(double length); 00531 00532 double x1;double y1;double z1; 00533 double x2;double y2;double z2; 00534 double x3;double y3;double z3; 00535 protected: 00536 void Stripeit(double length); 00537 void Newtriangles(double length); 00538 bool stripeit; 00539 double xa;double ya; 00540 double xb;double yb; 00541 }; 00542 00543 00544 class D3rectangle:public D3element{ 00545 public: 00546 virtual void InsertTGeoVolume(TGeoVolume *top,TGeoMedium *matVak,TGeoMedium *mat,TGeoManager *geom); 00547 virtual bool IntersectWithRay(double x0,double y0,double z0,double xdir,double ydir,double zdir,double &mul); 00548 virtual void SetNormTowards(double x,double y,double z,bool towards); 00549 virtual void GetReferencePoint(double &x,double &y, double &z); 00550 virtual double GetArea(); 00551 virtual void GetCenter(double &x,double &y,double &z); 00552 virtual double GetSelfPotential(); 00553 virtual inline double GetPotentialAt(double x,double y,double z); 00554 virtual double GetSelfDoubleLayerPotential(); 00555 virtual inline double GetDoubleLayerPotentialAt(double x,double y,double z); 00556 virtual void GetPotentialAndFieldAt(double x,double y,double z,double &pot,double &ex,double &ey,double &ez); 00557 virtual void GetRectangle(double *A,double *B,double *C,double *D,double *COL); 00558 virtual void rotate(double phi_,double theta_,double psi_); 00559 virtual void shift(double xs,double ys,double zs); 00581 D3rectangle(double xa_,double ya_,double xc_,double yc_, 00582 double x,double y,double z, 00583 double phi_,double theta_,double psi_,bool inversenorm=false,PD3element parent=NULL,bool subdivideInSqares=false); 00607 D3rectangle(double xa_,double ya_,double xb_,double yb_,double xc_,double yc_, 00608 double x,double y,double z, 00609 double phi_,double theta_,double psi_,bool inversenorm=false,PD3element parent=NULL); 00618 00619 D3rectangle(char *str,double x1_,double y1_,double z1_,double x2_,double y2_,double z2_,double x3_,double y3_,double z3_,double x4_,double y4_,double z4_, 00620 bool inversenorm=false,PD3element parent=NULL); 00621 ~D3rectangle(){}; 00622 virtual void createNewSubelements(double length); 00623 double x1;double y1;double z1; 00624 double x2;double y2;double z2; 00625 double x3;double y3;double z3; 00626 double x4;double y4;double z4; 00627 protected: 00628 double xa;double ya; 00629 double xb;double yb; 00630 double xc;double yc; 00631 }; 00632 00633 00634 class D3box:public D3element{ 00635 public: 00663 D3box(double xl_,double yl_,double zl_, 00664 double x,double y,double z, 00665 double phi_,double theta_,double psi_,bool inversenorm=false,PD3element parent=NULL): 00666 xl(xl_),yl(yl_),zl(zl_), 00667 D3element(x,y,z,phi_,theta_,psi_,inversenorm,false,parent){ 00668 subelement.push_back(new D3rectangle(xl,0,0,yl,0,0,0,0,0,0,!inversenorm,this));//triangle at the left bottom corner 00669 subelement.push_back(new D3rectangle(xl,0,0,yl,0,0,zl,0,0,0,inversenorm,this));//triangle at the right bottom corner 00670 subelement.push_back(new D3rectangle(xl,0,0,zl,0,0,0,0,90,0,inversenorm,this));//triangle at the right bottom corner 00671 subelement.push_back(new D3rectangle(xl,0,0,zl,0,yl,0,0,90,0,!inversenorm,this));//triangle at the right bottom corner 00672 subelement.push_back(new D3rectangle(yl,0,0,zl,0,0,0,0,90,90,!inversenorm,this));//triangle at the right bottom corner 00673 subelement.push_back(new D3rectangle(yl,0,0,zl,xl,0,0,0,90,90,inversenorm,this));//triangle at the right bottom corner 00674 00675 }; 00676 ~D3box(){}; 00677 00678 protected: 00679 double xl;double yl;double zl; 00680 }; 00681 00682 00683 class D3icosahedron:public D3element{ 00684 public: 00685 D3icosahedron(double edgelength,double refinement,double x,double y,double z,double phi_,double theta_,double psi_,bool inversenorm=false,PD3element parent=NULL):r(r),refinement(refinement), 00686 D3element(x,y,z,phi,theta,psi,inversenorm,false,parent){ 00687 double e=edgelength/2.,t=(1.+sqrt(5.))/2.*edgelength/2.; 00688 double VV[12][3]={{0,e,t},{e,t,0},{t,0,e}, 00689 {0,e,-t},{e,-t,0},{t,0,-e}, 00690 {0,-e,t},{-e,t,0},{-t,0,e}, 00691 {0,-e,-t},{-e,-t,0},{-t,0,-e}}; 00692 V=&VV; 00693 00694 00695 triangle( 0, 1, 2, true); 00696 triangle(1,2,5); 00697 triangle(4,6,2,true); 00698 00699 triangle(2,5,4,true); 00700 triangle(0,2,6,true); 00701 triangle(4,5,9,true); 00702 00703 triangle(4,9,10,true); 00704 triangle(4,6,10);//?? 00705 triangle(10,6,8); 00706 00707 triangle(8,6,0); 00708 triangle(8,7,0,true); 00709 triangle(0,1,7); 00710 00711 triangle(1,3,5,true); 00712 triangle(3,5,9); 00713 triangle(11,10,9,true); 00714 00715 triangle(11,10,8); 00716 triangle(11,8,7); 00717 triangle(11,7,3); 00718 triangle(11,3,9); 00719 00720 triangle(1 ,3,7); 00721 00722 00723 };//x --> eventuell fehler drin??, änderung von edgelength und refinement ändert nichts 00724 ~D3icosahedron(){}; 00725 void triangle(int a,int b,int c,bool inversenorm=false){ 00726 subelement.push_back(new D3triangle((*V)[a][0],(*V)[a][1],(*V)[a][2], 00727 (*V)[b][0],(*V)[b][1],(*V)[b][2], 00728 (*V)[c][0],(*V)[c][1],(*V)[c][2], 00729 inversenorm,this,false)); 00730 } 00731 protected: 00732 double r; 00733 double refinement; 00734 double (*V)[12][3]; 00735 }; 00736 00737 00738 class D3sphere:public D3element{ 00739 public: 00758 D3sphere(double r,double refinement,double x,double y,double z,double phi_,double theta_,double psi_,bool inversenorm=false,PD3element parent=NULL):r(r),refinement(refinement), 00759 D3element(x,y,z,phi,theta,psi,inversenorm,false,parent){ 00760 double edgelength=4./sqrt(10.+2.*sqrt(5.));//radius 1 00761 D3icosahedron *ico=new D3icosahedron(edgelength,0,0,0,0,0,0,inversenorm,NULL); 00762 ico->refine(refinement/r); 00763 int n=ico->GetAmountOfSubelements(); 00764 PD3element *el=new PD3element[n]; 00765 int cnt=0; 00766 ico->GetListOfBaseElements(el,cnt); 00767 00768 double a[3],b[3],c[3],anorm,bnorm,cnorm,color[3]; 00769 int col; 00770 for(col=0;col<n;col++){//geht alle Flaechenelemente durch 00771 if(dynamic_cast<D3triangle *>(el[col])) { 00772 el[col]->GetTriangle(a,b,c,color); 00773 anorm=sqrt(sqr(a[0])+sqr(a[1])+sqr(a[2])); 00774 bnorm=sqrt(sqr(b[0])+sqr(b[1])+sqr(b[2])); 00775 cnorm=sqrt(sqr(c[0])+sqr(c[1])+sqr(c[2])); 00776 subelement.push_back(new D3triangle(x+a[0]/anorm*r,y+a[1]/anorm*r,z+a[2]/anorm*r, 00777 x+b[0]/bnorm*r,y+b[1]/bnorm*r,z+b[2]/bnorm*r, 00778 x+c[0]/cnorm*r,y+c[1]/cnorm*r,z+c[2]/cnorm*r, 00779 inversenorm,this,false)); 00780 } 00781 } 00782 delete[] el; 00783 delete ico; 00784 00785 }; 00786 00787 00788 ~D3sphere(){}; 00789 protected: 00790 double r; 00791 double refinement; 00792 }; 00793 00794 class D3disk:public D3element{ 00795 public: 00815 00816 D3disk(double r_, 00817 double x,double y,double z, 00818 double phi_,double theta_,double psi_,int sectors_=32,bool subdivide=true,bool inversenorm=false,PD3element parent=NULL): 00819 r(r_),sectors(sectors_), 00820 D3element(x,y,z,phi_,theta_,psi_,inversenorm,false,parent){ 00821 int i; 00822 double step=2*pi/double(sectors); 00823 double dr=r*step; 00824 int subdivr=(r/dr); 00825 if(subdivr<1||!subdivide){ 00826 for(i=0;i<sectors;i++){ 00827 subelement.push_back(new D3triangle(r*cos(step*double(i)),r*sin(step*double(i)),r*cos(step*double(i+1)),r*sin(step*double(i+1)),0,0,0,0,0,0,inversenorm,this));//triangle at the left bottom corner 00828 } 00829 } 00830 else{ 00831 for(i=0;i<sectors;i++){ 00832 subelement.push_back(new D3triangle(r/double(subdivr)*cos(step*double(i)),r/double(subdivr)*sin(step*double(i)),r/double(subdivr)*cos(step*double(i+1)),r/double(subdivr)*sin(step*double(i+1)),0,0,0,0,0,0,inversenorm,this));//triangle at the left bottom corner 00833 int ring; 00834 for(ring=1;ring<subdivr;ring++){ 00835 double xa=r/double(subdivr)*double(ring)*cos(step*double(i)); 00836 double ya=r/double(subdivr)*double(ring)*sin(step*double(i)); 00837 double xd=r/double(subdivr)*double(ring)*cos(step*double(i+1)); 00838 double yd=r/double(subdivr)*double(ring)*sin(step*double(i+1)); 00839 double xb=r/double(subdivr)*double(ring+1)*cos(step*double(i)); 00840 double yb=r/double(subdivr)*double(ring+1)*sin(step*double(i)); 00841 double xc=r/double(subdivr)*double(ring+1)*cos(step*double(i+1)); 00842 double yc=r/double(subdivr)*double(ring+1)*sin(step*double(i+1)); 00843 subelement.push_back(new D3triangle(xb-xa,yb-ya,xc-xa,yc-ya,xa,ya,0,0,0,0,inversenorm,this)); 00844 subelement.push_back(new D3triangle(xc-xa,yc-ya,xd-xa,yd-ya,xa,ya,0,0,0,0,inversenorm,this)); 00845 } 00846 } 00847 00848 00849 } 00850 }; 00851 00852 ~D3disk(){}; 00853 protected: 00854 double r; 00855 double sectors; 00856 00857 }; 00858 00859 class D3cylinder:public D3element{ 00860 public: 00880 D3cylinder(double r_,double l_, 00881 double x,double y,double z, 00882 double phi_,double theta_,double psi_,int sectors_=12,bool subdivide=true,bool inversenorm=false,PD3element parent=NULL): 00883 r(r_),l(l_),sectors(sectors_), 00884 D3element(x,y,z,phi_,theta_,psi_,inversenorm,false,parent){ 00885 int i; 00886 double step=2*pi/double(sectors); 00887 double dr=r*sqrt(2.*(1.-cos(step))); 00888 int subdivr=(r/dr); 00889 if(subdivr<=1||!subdivide){ 00890 for(i=0;i<sectors;i++){ 00891 subelement.push_back(new D3triangle(r*cos(step*double(i)),r*sin(step*double(i)),r*cos(step*double(i+1)),r*sin(step*double(i+1)),0,0,0,0,0,0,!inversenorm,this));//triangle at the left bottom corner 00892 subelement.push_back(new D3triangle(r*cos(step*double(i)),r*sin(step*double(i)),r*cos(step*double(i+1)),r*sin(step*double(i+1)),0,0,l,0,0,0,inversenorm,this));//triangle at the left bottom corner 00893 double xb=r*cos(step*double(i)); 00894 double yb=r*sin(step*double(i)); 00895 subelement.push_back(new D3rectangle(dr,0,0,l,xb,yb,0,0,90,90.+180./pi*step*(double(i)+0.5),inversenorm,this,true)); 00896 00897 } 00898 00899 } 00900 else{ 00901 for(i=0;i<sectors;i++){ 00902 //center circle of front faces 00903 subelement.push_back(new D3triangle(r/double(subdivr)*cos(step*double(i)),r/double(subdivr)*sin(step*double(i)),r/double(subdivr)*cos(step*double(i+1)),r/double(subdivr)*sin(step*double(i+1)),0,0,0,0,0,0,!inversenorm,this));//triangle at the left bottom corner 00904 subelement.push_back(new D3triangle(r/double(subdivr)*cos(step*double(i)),r/double(subdivr)*sin(step*double(i)),r/double(subdivr)*cos(step*double(i+1)),r/double(subdivr)*sin(step*double(i+1)),0,0,l,0,0,0,inversenorm,this));//triangle at the left bottom corner 00905 int ring; 00906 double xb,yb; 00907 for(ring=1;ring<subdivr;ring++){ 00908 double xa=r/double(subdivr)*double(ring)*cos(step*double(i)); 00909 double ya=r/double(subdivr)*double(ring)*sin(step*double(i)); 00910 double xd=r/double(subdivr)*double(ring)*cos(step*double(i+1)); 00911 double yd=r/double(subdivr)*double(ring)*sin(step*double(i+1)); 00912 xb=r/double(subdivr)*double(ring+1)*cos(step*double(i)); 00913 yb=r/double(subdivr)*double(ring+1)*sin(step*double(i)); 00914 double xc=r/double(subdivr)*double(ring+1)*cos(step*double(i+1)); 00915 double yc=r/double(subdivr)*double(ring+1)*sin(step*double(i+1)); 00916 00917 //sectors of front faces as rectangles 00918 subelement.push_back(new D3rectangle(xb-xa,yb-ya,xc-xa,yc-ya,xd-xa,yd-ya,xa,ya,0,0,0,0,!inversenorm,this)); 00919 subelement.push_back(new D3rectangle(xb-xa,yb-ya,xc-xa,yc-ya,xd-xa,yd-ya,xa,ya,l,0,0,0,inversenorm,this)); 00920 //sectors of front faces as triangles (are documented out) 00921 // subelement.push_back(new D3triangle(xb-xa,yb-ya,xc-xa,yc-ya,xa,ya,0,0,0,0,!inversenorm,this)); 00922 // subelement.push_back(new D3triangle(xc-xa,yc-ya,xd-xa,yd-ya,xa,ya,0,0,0,0,!inversenorm,this)); 00923 00924 // subelement.push_back(new D3triangle(xb-xa,yb-ya,xc-xa,yc-ya,xa,ya,l,0,0,0,inversenorm,this)); 00925 // subelement.push_back(new D3triangle(xc-xa,yc-ya,xd-xa,yd-ya,xa,ya,l,0,0,0,inversenorm,this)); 00926 } 00927 //outer side element 00928 subelement.push_back(new D3rectangle(dr,0,0,l,xb,yb,0,0,90,90.+180./pi*step*(double(i)+0.5),inversenorm,this,true)); 00929 } 00930 00931 00932 } 00933 }; 00934 00935 ~D3cylinder(){}; 00936 protected: 00937 double r; 00938 double l; 00939 double sectors; 00940 00941 }; 00942 00943 00944 class D3tube:public D3element{ 00945 public: 00967 D3tube(double r_,double r2_,double l_, 00968 double x,double y,double z, 00969 double phi_,double theta_,double psi_,int sectors_=12,bool subdivide=true,bool inversenorm=false,PD3element parent=NULL): 00970 r(r_),r2(r2_),l(l_),sectors(sectors_), 00971 D3element(x,y,z,phi_,theta_,psi_,inversenorm,false,parent){ 00972 if(!subdivide){ 00973 cout <<"subdivide false not yet implemented"<<endl; 00974 } 00975 subdivide=true; 00976 int i; 00977 double step=2*pi/double(sectors); 00978 double dr=r*sqrt(2.*(1.-cos(step))); 00979 int subdivr=((r-r2)/dr); 00980 if(subdivr<=1) subdivr=2; 00981 if(subdivr<=1||!subdivide){ 00982 for(i=0;i<sectors;i++){ 00983 subelement.push_back(new D3triangle(r*cos(step*double(i)),r*sin(step*double(i)),r*cos(step*double(i+1)),r*sin(step*double(i+1)),0,0,0,0,0,0,!inversenorm,this));//triangle at the left bottom corner 00984 subelement.push_back(new D3triangle(r*cos(step*double(i)),r*sin(step*double(i)),r*cos(step*double(i+1)),r*sin(step*double(i+1)),0,0,l,0,0,0,inversenorm,this));//triangle at the left bottom corner 00985 double xb=r*cos(step*double(i)); 00986 double yb=r*sin(step*double(i)); 00987 subelement.push_back(new D3rectangle(dr,0,0,l,xb,yb,0,0,90,90.+180./pi*step*(double(i)+0.5),inversenorm,this,true)); 00988 00989 } 00990 00991 } 00992 else{ 00993 for(i=0;i<sectors;i++){ 00994 //center circle of front faces 00995 //subelement.push_back(new D3triangle(r/double(subdivr)*cos(step*double(i)),r/double(subdivr)*sin(step*double(i)),r/double(subdivr)*cos(step*double(i+1)),r/double(subdivr)*sin(step*double(i+1)),0,0,0,0,0,0,!inversenorm,this));//triangle at the left bottom corner 00996 //subelement.push_back(new D3triangle(r/double(subdivr)*cos(step*double(i)),r/double(subdivr)*sin(step*double(i)),r/double(subdivr)*cos(step*double(i+1)),r/double(subdivr)*sin(step*double(i+1)),0,0,l,0,0,0,inversenorm,this));//triangle at the left bottom corner 00997 int ring; 00998 double xb,yb; 00999 for(ring=0;ring<subdivr;ring++){ 01000 double rr=r2+(r-r2)/double(subdivr)*double(ring); 01001 double rr2=r2+(r-r2)/double(subdivr)*double(ring+1); 01002 double xa=rr*cos(step*double(i)); 01003 double ya=rr*sin(step*double(i)); 01004 double xd=rr*cos(step*double(i+1)); 01005 double yd=rr*sin(step*double(i+1)); 01006 xb=rr2*cos(step*double(i)); 01007 yb=rr2*sin(step*double(i)); 01008 double xc=rr2*cos(step*double(i+1)); 01009 double yc=rr2*sin(step*double(i+1)); 01010 01011 //sectors of front faces as rectangles 01012 subelement.push_back(new D3rectangle(xb-xa,yb-ya,xc-xa,yc-ya,xd-xa,yd-ya,xa,ya,0,0,0,0,!inversenorm,this)); 01013 subelement.push_back(new D3rectangle(xb-xa,yb-ya,xc-xa,yc-ya,xd-xa,yd-ya,xa,ya,l,0,0,0,inversenorm,this)); 01014 //sectors of front faces as triangles (are documented out) 01015 // subelement.push_back(new D3triangle(xb-xa,yb-ya,xc-xa,yc-ya,xa,ya,0,0,0,0,!inversenorm,this)); 01016 // subelement.push_back(new D3triangle(xc-xa,yc-ya,xd-xa,yd-ya,xa,ya,0,0,0,0,!inversenorm,this)); 01017 01018 // subelement.push_back(new D3triangle(xb-xa,yb-ya,xc-xa,yc-ya,xa,ya,l,0,0,0,inversenorm,this)); 01019 // subelement.push_back(new D3triangle(xc-xa,yc-ya,xd-xa,yd-ya,xa,ya,l,0,0,0,inversenorm,this)); 01020 } 01021 //outer side element 01022 subelement.push_back(new D3rectangle(dr,0,0,l,xb,yb,0,0,90,90.+180./pi*step*(double(i)+0.5),inversenorm,this,true)); 01023 //inner side element 01024 01025 xb=r2*cos(step*double(i)); 01026 yb=r2*sin(step*double(i)); 01027 double xb2=r2*cos(step*double(i+1)); 01028 double yb2=r2*sin(step*double(i+1)); 01029 double dr2=sqrt((xb-xb2)*(xb-xb2)+(yb-yb2)*(yb-yb2)); 01030 subelement.push_back(new D3rectangle(dr2,0,0,l,xb,yb,0,0,90,90.+180./pi*step*(double(i)+0.5),!inversenorm,this,true)); 01031 } 01032 01033 01034 } 01035 }; 01036 ~D3tube(){}; 01037 protected: 01038 double r; 01039 double r2; 01040 double l; 01041 double sectors; 01042 01043 }; 01044 01045 01047 // Electron charge 01048 #define eee 1.602e-19 01049 //Mass of Calcium 01050 #define mmm (40.078*1.66e-27) 01051 01052 class D3world:public D3element{ 01053 01054 public: 01055 bool rangeerror; 01056 D3world(unsigned long num1,unsigned long num2,unsigned long num3,unsigned long num4,unsigned long d,unsigned long m,unsigned long y); //x 01061 01062 D3world(const char *_cachefilename,double tol=0.001,int maxit=32,int numMom=2,int numLev=4,double spaceunit=0.001,int segmentation=1000); 01063 ~D3world(); 01064 void insert(D3electrode *el);//x 01065 void GetListOfBaseElements(PD3element *el,int *electrodeIndexLimit,int &cnt); 01066 void SymmetrizeCharges(int axis=1,double epsilon=0.00001,bool ignoremirror=false);//x=1;y=2;z=3 call it after solve!//x 01067 void solve();//x ?? 01071 double calc(double x,double y,double z);//x 01072 void calc(double xxx,double y,double z,double &pot,double &feldx,double &feldy,double &feldz);//x 01073 01074 void calc(double xmin,double xmax,int nx,double ymin,double ymax,int ny,double zmin,double zmax,int nz); 01075 void calc(int xxxnum,double *xxx,double *Potential); 01076 void calc(int xxxnum,double *xxx,double *Potential,double *FieldX,double *FieldY,double *FieldZ); 01077 double calc_slow(double x,double y,double z); //x 01078 void calc_slow(double xmin,double xmax,int nx,double ymin,double ymax,int ny,double zmin,double zmax,int nz); //x 01079 void calc_slow(int xxxnum,double *xxx,double *Potential,double *FieldX,double *FieldY,double *FieldZ); 01080 void SetScalePostCalc(double xscale2,double yscale2,double zscale2); 01081 bool IsEqualSurfaceElement(int amountOfVertices,double eps,double *x,double *X); 01082 void draw(); //x 01083 void Dcentroid(int shape,double * pc,double* xcout); 01084 void save(char *fname); 01085 bool load(char *fname);//true if saved data is compatible 01086 int cut(int n,int max); 01087 void savecalc(char *fname); 01088 bool loadcalc(char *fname);//true if saved data is compatible 01089 void AssignColors(); 01090 unsigned long update_adler32double(unsigned long old, double *buf, unsigned long len,int ignorebytes=3); 01091 void RefreshChecksum(); 01092 void exportGeometry(const char *fname);//x 01098 void propagateForwardVerlet(double x[3],double v[3],double h,double qDivM=eee/mmm,bool onedim=false);//x 01099 void propagateForwardEuler(double x[3],double v[3],double h,double qDivM=eee/mmm);//x 01100 void propagateForwardVerletRotSymX(double x[3],double v[3],double h,double qDivM=eee/mmm);//x 01101 void propagateForwardVerletRotSymY(double x[3],double v[3],double h,double qDivM=eee/mmm);//x 01102 void propagateForwardVerletRotSymZ(double x[3],double v[3],double h,double qDivM=eee/mmm);//x 01103 void calc2(int xxxnum,double *xxx,double *Potential,double *FieldX,double *FieldY,double *FieldZ); 01104 void calc_slow2(int xxxnum,double *xxx,double *Potential,double *FieldX,double *FieldY,double *FieldZ); 01105 01106 protected: 01107 double xscale; 01108 double yscale; 01109 double zscale; 01110 int totalrot; 01111 void exch(double &a,double &b); 01112 void RotMirrorSurfaceElement(int axis,int sym,double *xx); 01113 bool refreshchecksum; 01114 int currentel; 01115 int segmentation; 01116 double spaceunit; 01117 char cachefilename[2000]; 01118 //caching stuff 01119 bool docache; 01120 double xmin; 01121 double xmax; 01122 int nx; 01123 double ymin; 01124 double ymax; 01125 int ny; 01126 double zmin; 01127 double zmax; 01128 int nz; 01129 double* potcache; 01130 double* feldxcache; 01131 double* feldycache; 01132 double* feldzcache; 01133 //end of caching stuff 01134 int size; 01135 int nlhs; 01136 int nrhs; 01137 int numMom; 01138 int numLev; 01139 int i; 01140 int j; 01141 char *shapechar; 01142 double *x; 01143 double *poten; 01144 double *dbydnpotenAll; 01145 double *dbydbpoten; 01146 double *xcoll; 01147 double *xnrm; 01148 double *lhsvect; 01149 double *rhsvect; 01150 int *shape; 01151 int *type; 01152 int *dtype; 01153 int *rhstype; 01154 int *lhstype; 01155 int *rhsindex; 01156 int *lhsindex; 01157 int job; 01158 int fljob; 01159 double error; 01160 double max_diri; 01161 double ave_diri; 01162 double max_neum; 01163 double ave_neum; 01164 double cnt_diri; 01165 double cnt_neum; 01166 /* Set the tolerance and max iterations for GMRES called by fastlap. */ 01167 double tol; 01168 int maxit; 01169 int numit; 01170 unsigned long checksum; 01171 unsigned long checksumcalc; 01172 01173 double *f; 01174 double *dfdnAll; 01175 int *electrodeIndexLimit; 01176 PD3element *el; 01177 D3electrode **electrodes; 01178 int n; 01179 int amountOfElectrodes; 01180 unsigned long update_adler32(unsigned long old, unsigned char *buf, unsigned long len); 01181 01182 }; 01184 class D3slowworld:public D3element{ 01185 public: 01186 D3slowworld(); 01187 ~D3slowworld(); 01188 void insert(D3electrode *el); 01189 void GetListOfBaseElements(PD3element *el,int *electrodeIndexLimit,int &cnt); 01190 void solve(); 01191 double calc(double x,double y,double z); 01192 void draw(); 01193 01194 protected: 01195 double *f; 01196 double *dfdnAll; 01197 int *electrodeIndexLimit; 01198 PD3element *el; 01199 int n; 01200 int amountOfElectrodes; 01201 }; 01202 01203 extern calcD3triangle calc; 01204 01205 01206 class D3thicktriangle:public D3element{ 01207 public: 01227 D3thicktriangle(double xa,double ya,double xb,double yb,double htriangle, 01228 double x,double y,double z, 01229 double phi,double theta,double psi,bool inversenorm=false,PD3element parent=NULL): 01230 xa(xa),ya(ya),xb(xb),yb(yb),htriangle(htriangle), 01231 D3element(x,y,z,phi,theta,psi,inversenorm,false,parent){ 01232 subelement.push_back(new D3triangle(xa,ya,xb,yb,0,0,0,0,0,0,!inversenorm,this)); 01233 subelement.push_back(new D3triangle(xa,ya,xb,yb,0,0,htriangle,0,0,0,inversenorm,this)); 01234 01235 double l1=sqrt(sqr(xa)+sqr(ya)); 01236 double w1=atan(ya/xa); 01237 subelement.push_back(new D3rectangle(l1,0,0,htriangle,0,0,0,0,90,w1/pi*180,!inversenorm,this)); 01238 01239 double l2=sqrt(sqr(xb)+sqr(yb)); 01240 double w2=atan(yb/xb); 01241 subelement.push_back(new D3rectangle(l2,0,0,htriangle,0,0,0,0,90,w2/pi*180,inversenorm,this)); 01242 01243 double l3=sqrt(sqr(xb-xa)+sqr(yb-ya)); 01244 double w3=atan((xb-xa)/(yb-ya)); 01245 01246 subelement.push_back(new D3rectangle(l3,0,0,htriangle,xa,ya,0,0,90,90-w3/pi*180,inversenorm,this)); 01247 01248 01249 01250 }; 01251 01252 ~D3thicktriangle(){}; 01253 protected: 01254 double xa;double ya;double xb;double yb;double htriangle; 01255 }; 01256 01257 01258 01259 #endif
