Метод магнитостатики конечных элементов на треугольниках, страница 5

ig=(long*)calloc(n+1,sizeof(long));

for(long i=0;i<=n;i++)

ig[i]=a.ig[i];

jg=(long*)calloc(ig[n]-1,sizeof(long));

gg=(double*)calloc(ig[n]-1,sizeof(double));

for(long i=0;i<ig[n]-1;i++)

{

jg[i]=a.jg[i];

gg[i]=a.gg[i];

}

}

//-------------------------------------------------------MATRIX::~MATRIX()

{

if(d!=NULL)free(d);

if(ig!=NULL)free(ig);

if(jg!=NULL)free(jg);

if(gg!=NULL)free(gg);

}

//---------------------------------------------------------void MATRIX::formier_profil(field my)

{

//формируем временный список для профиля

//В нём будут хранится какие эл-ты ненукевые

//Список temp[i] говорит о наличие эл-тов (i,el списка)

//Таким образом, кол-во эл-тов в temp[i] - данные массива ig

//Сами эл-ты списка temp[i] - данные массива jg

list_long *temp=(list_long*)calloc(my.m_knot.n,sizeof(list_long));

for(long i=0;i<my.m_loc.n;i++)

{

temp[my.m_loc.mas[i].num_knot[2]-1].add_el(my.m_loc.mas[i].num_knot[0]);

temp[my.m_loc.mas[i].num_knot[2]-1].add_el(my.m_loc.mas[i].num_knot[1]);

temp[my.m_loc.mas[i].num_knot[1]-1].add_el(my.m_loc.mas[i].num_knot[0]);

}

//Определение параметров для матрицы

n=my.m_knot.n;

d=(double*)calloc(n,sizeof(double));

//for(long i=0;i<n;i++) d[i]=0;

ig=(long*)calloc(n+1,sizeof(long));

//Заполнение массива ig

ig[0]=1;

long amm=0,tmp=0;

for(long i=1;i<=n;i++)

{

tmp=temp[i-1].get_ammount();

amm+=tmp;

ig[i]=ig[i-1]+tmp;

}

jg=(long*)calloc(amm,sizeof(long));

gg=(double*)calloc(amm,sizeof(double));

//Заполнение массива jg

for(long i=0,j=0;i<amm&&j<my.m_knot.n;)

{

temp[j].copy_to_mas(&jg[i]);

//printf("i=%d j=%d\n",i,j);

i+=temp[j].get_ammount();

j++;

}

//Освобождение временной памяти

for(long i=0;i<my.m_knot.n;i++)

temp[i].free_list();

//Профиль готов!!!

}

//----------------------------------------------------------------void MATRIX::add(long i,long j,double x)

{

if(i<j)

{long temp=i;i=j;j=temp;}

long k;

int flag=1;

for(k=ig[i]-1;k<ig[i+1]-1&&flag;k++)

if(jg[k]==j+1) flag=0;

if(!flag)

gg[k-1]+=x;

else

printf("ERR in MATRIX::add\n");

}

//-------------------------------------------------------void SLAU::add_loc_in_global(LOCAL_MATRIX my,VECTOR f_loc,loc lc)

{

long l,k;

for(long i=0;i<M;i++)

{

l=lc.num_knot[i]-1;

for(long j=i+1;j<M;j++)

{

k=lc.num_knot[j]-1;

add(l,k,my.mas[i][j]);

}

f.vect[l]+=f_loc.vect[i];

d[l]+=my.mas[i][i];

}

}

//--------------------------------------------------------void SLAU::formier_slau(field my)

{

formier_profil(my);

x.call_memory(n);

f.call_memory(n);

VECTOR f_loc(M),ff(M);//локальный вектор f

LOCAL_MATRIX C,loc_c,loc_b,loc;

C.init_c();// шаблон локальной матрицы

//пробегаем по локальным эл-там

for(long k=0;k<my.m_loc.n;k++)

{

ff.ins_el_all(my.get_toku(k));

loc_c=my.mes(k)*C; //формирование локальной матрицы массы

f_loc=loc_c*ff; //Формирование локального вектора правой части

loc_b.init_b_loc(my,k);

loc=loc_b;

add_loc_in_global(loc,f_loc,my.m_loc.mas[k]); //добавление в глобал. матрицу*/

}

//Учёт краевых условий

long tmp;

for(long i=0;i<my.m_kt1.n;i++)

{

tmp=my.m_kt1.mas[i].num-1;

d[tmp]=MAX;

f.ins_el(my.m_kt1.mas[i].a*MAX,tmp);

}

/*VECTOR xx(n);

double aa;

int des=open("d:\\file\\v2.dat",O_RDONLY|O_BINARY);

for(long i=0;i<n;i++)

{

read(des,&aa,sizeof(double));

xx.ins_el(aa,i);

}

close(des);

VECTOR f2=*this*xx;

double ss=(f2-f).norm()/f.norm();

if(ss<EPS) {printf("MATRIX OK\n");getch();}*/

}

//--------------------------------------------------------------------------void SLAU::MSG()

{

MATRIX s;

SST(&s);

VECTOR r(n),z(n),h(n),h2(n),h3(n);

x.ins_el_all(0);

r=f-*this*x;

z=s.solution_x(r);

double f_norm=f.norm();

double tmp;

double e2=r.norm()/f_norm;

double a,b;

long k;

h2=s.solution_x(r); //(SST)^-1*r в h2

tmp=h2*r;

for(k=1;e2>EPS;k++)

{

//a

h=*this*z;//A*z

a=tmp/(h*z);

//r

r=r-a*h;

//x

x=x+a*z;

//b

b=1./tmp;

h2=solution_x(r); //(SST)^-1*r в h2

tmp=h2*r;

b*=tmp;

//z

z=h2+b*z;

e2=r.norm()/f_norm;

double eee=fabs(e2-(f-*this*x).norm()/f.norm());

if(eee>EPS)

{printf("ERROR eps-e\n");getch();}