assemble_matrix_A.h

Go to the documentation of this file.

/*
   This file belongs to Aeneas. Aeneas is a GNU package released under GPL 3 license.
   This code is a simulator for Submicron 3D Semiconductor Devices. 
   It implements the Monte Carlo transport in 3D tetrahedra meshes
   for the simulation of the semiclassical Boltzmann equation for both electrons.
   It also includes all the relevant quantum effects for nanodevices.

   Copyright (C) 2007 Jean Michel Sellier <sellier@dmi.unict.it>
 
   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 3, or (at your option)
   any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program. If not, see <http://www.gnu.org/licenses/>.
*/

// Created on : 06 june 2007, Siracusa, Jean Michel Sellier
// Last modified : 16 august 2007, Siracusa, Jean Michel Sellier

// Definition of the A matrix in a row-indexed format for the resolution 
// of the discretised Poisson equation by means of finite
// element. The boundary conditions for the potential are of the type
// Dirichlet/Neumann mixed conditions. In order to take into account the
// non-homogeneous Dirichlet boundary conditions we implement the
// well-known Iron technique. For more informations see french book
// pag.336 par.9.5.5.

// This routine automaticaly sets the field continuity on the heterojunctions

void assemble_matrix_A(void)
{
// int dum;
 int m,l,k1,k2,ni,nj;
 double xii[3][3]; // i=1..dim, j=1..lg
 double *s[3][3]; // i=1..dim, j=1..dim
 double sigma[3][3]; // i=1..dim, j=1..dim
 double *A;

 int i,j;
// int row_start,row_end;

 for(i=0;i<3;i++) for(j=0;j<3;j++){xii[i][j]=0.; sigma[i][j]=0.;}
// memset(&xii,0,sizeof(xii));
// memset(&sigma,0,sizeof(sigma));

// set the sigma matrix
 for(i=0;i<3;i++) for(j=0;j<3;j++) s[i][j]=malloc((Ne+1)*sizeof(double));
 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...
 
 for(i=0;i<Ne;i++){
   if(i_dom[i]!=SIO2) s[0][0][i]=(epslon0*EPSR[i_dom[i]]);
   if(i_dom[i]==SIO2) s[0][0][i]=(epslon0*EPSR_SiO2);
   
   s[0][1][i]=0.;
   s[0][2][i]=0.;

   s[1][0][i]=0.;
   s[1][1][i]=s[0][0][i];
   s[1][2][i]=0.;

   s[2][0][i]=0.;
   s[2][1][i]=0.;
   s[2][2][i]=s[0][0][i];
 }
 // create the matrix
 printf("\nAllocating memory for the sparse Poisson matrix...\n");
 
 A = malloc((Ng+1)*(Ng+1)*sizeof(double)); //A[1..Ng][1..Ng];
 if(A==NULL){
   printf("assemble_matrix_A error : not enought memory for A allocation...\n");
   exit(0);       
 }
 printf("The Poisson matrix has been allocated with success.\n");
 printf("Matrix global sizes: sizeof(A[%d][%d]) = %d bytes\n\n", 
         Ng+1, Ng+1, (Ng+1)*(Ng+1)*sizeof(double));
// REMEMBER !!! A[i][j]=A[i*ncol+j];


// We set the elements of the matrix A[i][j] for i=1..Ng, j=1..Ng
// For more informations see the french book, pag.325, par.9.3.1 (cas 1)
 printf("setting the values for the Poisson matrix...\n");
 for(m=0;m<Ne;m++){
  double x1;
//  double x2=0.0;
//  memset(&xii,0.0,sizeof(xii));
  for(l=0;l<lg;l++) for(k1=0;k1<dim;k1++) xii[k1][l]=0.;
  for(l=0;l<lg;l++){
   for(k1=0;k1<dim;k1++){
// Evaluer les coordonees cartesiennes du point de Gauss \xi_l
     {int n; for(n=0;n<ng;n++) xii[k1][l]+=coord[k1][noeud_geo[n][m]-1]*base_geo[n][l];}
// Evaluer la valeur de sigma=1./(epslon0*epslon_Si) * I sur les noeuds
     for(k2=0;k2<dim;k2++){
       int n; for(n=0;n<nf;n++) sigma[k1][k2]+=s[k1][k2][m]*base_ref[n][l];
     }
   }
   for(ni=0;ni<nf;ni++){
    i=noeud_geo[ni][m]-1;
    for(nj=0;nj<nf;nj++){
     j=noeud_geo[nj][m]-1;
     {int k1i,k2i; for(k1i=0;k1i<dim;k1i++)
      for(k2i=0;k2i<dim;k2i++)
       x1+=d_base_K[k1i][nj][m]*s[k1i][k2i][m]*d_base_K[k2i][ni][m]; // Laplacien
      x1*=poids_K[l][m];
//      memset(&sigma,0.0,sizeof(sigma));
      {int hi,hj; for(hi=0;hi<3;hi++) for(hj=0;hj<3;hj++) sigma[hi][hj]=0.;}
     }
     A[i*Ng+j]+=x1;
    x1=0.0;
    }
   }
  }
 }
// non-homogeneous Dirichlet boundary conditions
// *****
 for(m=0;m<Ne;m++)
  for(ni=0;ni<nf;ni++){
    double x2;
    i=noeud_geo[ni][m]-1;
    if(i_front[i]==OHMIC || i_front[i]==SCHOTTKY){
     x2=1./EPS;
     A[i*Ng+i] += x2;
     x2=0.0;
    }
  }
  printf("setting done.\n\n");
// *****

 printf("Converting full storage mode into row-indexed sparse storage mode...\n");
#include "../linear_solver/sprsin.h"
 printf("Sparse storage mode done...\n\n");

// Deallocation of matrix A
 free(A);
 for(i=0;i<3;i++) for(j=0;j<3;j++)  free(s[i][j]);
 printf("Poisson matrix deallocated with success...\n\n");
}

---

© sourcejam.com 2005-2008