Main Page | Directories | File List | File Members

assemble_matrix_A.h File Reference

#include "../linear_solver/sprsin.h"

Go to the source code of this file.

Functions

void assemble_matrix_A (void)


Function Documentation

void assemble_matrix_A void   ) 
 

Definition at line 37 of file assemble_matrix_A.h.

References base_geo, base_ref, coord, d_base_K, dim, EPS, epslon0, EPSR, EPSR_SiO2, free(), i_dom, i_front, lg, Ne, nf, ng, Ng, noeud_geo, OHMIC, poids_K, printf(), SCHOTTKY, and SIO2.

Referenced by main().

00038 {
00039 // int dum;
00040  int m,l,k1,k2,ni,nj;
00041  double xii[3][3]; // i=1..dim, j=1..lg
00042  double *s[3][3]; // i=1..dim, j=1..dim
00043  double sigma[3][3]; // i=1..dim, j=1..dim
00044  double *A;
00045 
00046  int i,j;
00047 // int row_start,row_end;
00048 
00049  for(i=0;i<3;i++) for(j=0;j<3;j++){xii[i][j]=0.; sigma[i][j]=0.;}
00050 // memset(&xii,0,sizeof(xii));
00051 // memset(&sigma,0,sizeof(sigma));
00052 
00053 // set the sigma matrix
00054  for(i=0;i<3;i++) for(j=0;j<3;j++) s[i][j]=malloc((Ne+1)*sizeof(double));
00055  for(i=0;i<3;i++) for(j=0;j<3;j++) for(l=0;l<Ne;l++) s[i][j][l]=0.; // just in case...
00056  
00057  for(i=0;i<Ne;i++){
00058    if(i_dom[i]!=SIO2) s[0][0][i]=(epslon0*EPSR[i_dom[i]]);
00059    if(i_dom[i]==SIO2) s[0][0][i]=(epslon0*EPSR_SiO2);
00060    
00061    s[0][1][i]=0.;
00062    s[0][2][i]=0.;
00063 
00064    s[1][0][i]=0.;
00065    s[1][1][i]=s[0][0][i];
00066    s[1][2][i]=0.;
00067 
00068    s[2][0][i]=0.;
00069    s[2][1][i]=0.;
00070    s[2][2][i]=s[0][0][i];
00071  }
00072  // create the matrix
00073  printf("\nAllocating memory for the sparse Poisson matrix...\n");
00074  
00075  A = malloc((Ng+1)*(Ng+1)*sizeof(double)); //A[1..Ng][1..Ng];
00076  if(A==NULL){
00077    printf("assemble_matrix_A error : not enought memory for A allocation...\n");
00078    exit(0);       
00079  }
00080  printf("The Poisson matrix has been allocated with success.\n");
00081  printf("Matrix global sizes: sizeof(A[%d][%d]) = %d bytes\n\n", 
00082          Ng+1, Ng+1, (Ng+1)*(Ng+1)*sizeof(double));
00083 // REMEMBER !!! A[i][j]=A[i*ncol+j];
00084 
00085 
00086 // We set the elements of the matrix A[i][j] for i=1..Ng, j=1..Ng
00087 // For more informations see the french book, pag.325, par.9.3.1 (cas 1)
00088  printf("setting the values for the Poisson matrix...\n");
00089  for(m=0;m<Ne;m++){
00090   double x1;
00091 //  double x2=0.0;
00092 //  memset(&xii,0.0,sizeof(xii));
00093   for(l=0;l<lg;l++) for(k1=0;k1<dim;k1++) xii[k1][l]=0.;
00094   for(l=0;l<lg;l++){
00095    for(k1=0;k1<dim;k1++){
00096 // Evaluer les coordonees cartesiennes du point de Gauss \xi_l
00097      {int n; for(n=0;n<ng;n++) xii[k1][l]+=coord[k1][noeud_geo[n][m]-1]*base_geo[n][l];}
00098 // Evaluer la valeur de sigma=1./(epslon0*epslon_Si) * I sur les noeuds
00099      for(k2=0;k2<dim;k2++){
00100        int n; for(n=0;n<nf;n++) sigma[k1][k2]+=s[k1][k2][m]*base_ref[n][l];
00101      }
00102    }
00103    for(ni=0;ni<nf;ni++){
00104     i=noeud_geo[ni][m]-1;
00105     for(nj=0;nj<nf;nj++){
00106      j=noeud_geo[nj][m]-1;
00107      {int k1i,k2i; for(k1i=0;k1i<dim;k1i++)
00108       for(k2i=0;k2i<dim;k2i++)
00109        x1+=d_base_K[k1i][nj][m]*s[k1i][k2i][m]*d_base_K[k2i][ni][m]; // Laplacien
00110       x1*=poids_K[l][m];
00111 //      memset(&sigma,0.0,sizeof(sigma));
00112       {int hi,hj; for(hi=0;hi<3;hi++) for(hj=0;hj<3;hj++) sigma[hi][hj]=0.;}
00113      }
00114      A[i*Ng+j]+=x1;
00115     x1=0.0;
00116     }
00117    }
00118   }
00119  }
00120 // non-homogeneous Dirichlet boundary conditions
00121 // *****
00122  for(m=0;m<Ne;m++)
00123   for(ni=0;ni<nf;ni++){
00124     double x2;
00125     i=noeud_geo[ni][m]-1;
00126     if(i_front[i]==OHMIC || i_front[i]==SCHOTTKY){
00127      x2=1./EPS;
00128      A[i*Ng+i] += x2;
00129      x2=0.0;
00130     }
00131   }
00132   printf("setting done.\n\n");
00133 // *****
00134 
00135  printf("Converting full storage mode into row-indexed sparse storage mode...\n");
00136 #include "../linear_solver/sprsin.h"
00137  printf("Sparse storage mode done...\n\n");
00138 
00139 // Deallocation of matrix A
00140  free(A);
00141  for(i=0;i<3;i++) for(j=0;j<3;j++)  free(s[i][j]);
00142  printf("Poisson matrix deallocated with success...\n\n");
00143 }


© sourcejam.com 2005-2008