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

double x3=m_knot.mas[m_loc.mas[i].num_knot[2]-1].x;

double y1=m_knot.mas[m_loc.mas[i].num_knot[0]-1].y;

double y2=m_knot.mas[m_loc.mas[i].num_knot[1]-1].y;

double y3=m_knot.mas[m_loc.mas[i].num_knot[2]-1].y;

return fabs((x2*y3+x1*y2+x3*y1)-(x2*y1+x1*y3+x3*y2))/2.;

}

//--------------------------------------------double field::get_mu(long i)

{

long num=m_loc.mas[i].num_mat;

for(long j=0;j<m_mat.n;j++)

if(num==m_mat.mas[j].num) return m_mat.mas[j].mu;

printf("ERR GET_MU\n");

return NULL;

}

//--------------------------------------------double field::get_toku(long i)

{

long num=m_loc.mas[i].num_mat;

for(long j=0;j<m_mat.n;j++)

if(num==m_mat.mas[j].num) return m_mat.mas[j].toku;

printf("ERR GET_TOKU\n");

return NULL;

}

//-------function of list_long---------------list_long::list_long()

{el=0;next=NULL;}

//----------------------------------void list_long::add_el(long ell)

{

list_long *q=next,*p;

el++;

if(next==NULL)

{

next=new list_long;

next->el=ell;

}

else

{

if(ell<next->el)

{

p=next;

next=new list_long;

next->el=ell;

next->next=p;

}

else

{

while(q!=NULL&&q->el<=ell)

{p=q; q=q->next;}

if(p->el!=ell)

{

p->next=new list_long;

p->next->el=ell;

p->next->next=q;

}

else el--;

}

}

}

//--------------------------------------void list_long::del_el(long ell)

{

list_long *q=next,*p;

el--;

if(q->el==ell)

{

next=next->next;

delete q;

}

else

{

while(q!=NULL&&q->el<ell)

{p=q;q=q->next;}

if(q!=NULL&&q->el==ell)

{

p->next=q->next;

delete q;

}

else el++;

}

}

//-----------------------------------void list_long::free_list()

{

list_long *q=next,*p;

while(q!=NULL)

{p=q;q=q->next;delete p;}

}

//----------------------------------void list_long::print_list()

{

list_long *q=next;

while(q!=NULL)

{printf("%5d ",q->el);q=q->next;}

printf("\tn=%d\n",get_ammount());

}

//---------------------------------------long list_long::get_ammount()

{return el;}

void list_long::copy_to_mas(long *mas)

{

list_long *q=next;

for(long i=0;q!=NULL;i++)

{

mas[i]=q->el;

q=q->next;

}

}

//-------------function of LOCAL_MATRIX-----------void LOCAL_MATRIX::init_c()

{

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

{

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

mas[i][j]=1./12;

mas[i][i]*=2;

}

}

//--------------------------------------------------------------void LOCAL_MATRIX::init_b_loc(field my,long i)

{

double koeff;

koeff=(my.mes(i)*4*my.get_mu(i)*mu0);

double x1=my.m_knot.mas[my.m_loc.mas[i].num_knot[0]-1].x;

double x2=my.m_knot.mas[my.m_loc.mas[i].num_knot[1]-1].x;

double x3=my.m_knot.mas[my.m_loc.mas[i].num_knot[2]-1].x;

double y1=my.m_knot.mas[my.m_loc.mas[i].num_knot[0]-1].y;

double y2=my.m_knot.mas[my.m_loc.mas[i].num_knot[1]-1].y;

double y3=my.m_knot.mas[my.m_loc.mas[i].num_knot[2]-1].y;

double l1[]={y2-y3,y3-y1,y1-y2};

double l2[]={x3-x2,x1-x3,x2-x1};

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

for(long k=0;k<=l;k++)

mas[l][k]=(l1[l]*l1[k]+l2[l]*l2[k])/koeff;

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

for(long k=l+1;k<M;k++)

mas[l][k]=mas[k][l];

}

//--------------------------------------------------------------LOCAL_MATRIX LOCAL_MATRIX::operator+(LOCAL_MATRIX ob2)

{

LOCAL_MATRIX temp;

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

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

temp.mas[i][j]=mas[i][j]+ob2.mas[i][j];

return temp;

}

//----------------------------------------------------------LOCAL_MATRIX operator*(double a,LOCAL_MATRIX ob)

{

LOCAL_MATRIX temp;

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

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

temp.mas[i][j]=a*ob.mas[i][j];

return temp;

}

//---------------------------------------------------------LOCAL_MATRIX LOCAL_MATRIX::operator=(LOCAL_MATRIX ob2)

{

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

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

mas[i][j]=ob2.mas[i][j];

return *this;

}

//---------------------------------------------------------VECTOR operator*(LOCAL_MATRIX m,VECTOR v)

{

if(v.n!=M) {printf("ERROR IN LOCAL_MATRIX*VECTOR\n");return NULL;}

VECTOR temp(M);

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

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

temp.vect[i]+=m.mas[i][j]*v.vect[j];

return temp;

}

//-------------function of VECTOR--------------------------VECTOR::VECTOR()