Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5820) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5821) @@ -1275,5 +1275,4 @@ int i; - bool found = false; Input *input = NULL; @@ -2395,46 +2394,27 @@ void Tria::CreateKMatrixBalancedthickness_CG(Mat Kgg){ - /* local declarations */ - int i,j; - - /* node data: */ + /*Constants*/ const int numdof=NDOF1*NUMVERTICES; - double xyz_list[NUMVERTICES][3]; - int* doflist=NULL; - - /* gaussian points: */ - int ig; - GaussTria *gauss=NULL; - - /* matrices: */ - double L[NUMVERTICES]; - double B[2][NUMVERTICES]; - double Bprime[2][NUMVERTICES]; - double DL[2][2]={0.0}; - double DLprime[2][2]={0.0}; - double DL_scalar; - double Ke_gg[numdof][numdof]={0.0};//local element stiffness matrix - double Ke_gg_gaussian[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point. - double Ke_gg_thickness1[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point. - double Ke_gg_thickness2[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point. - - double Jdettria; - - /*input parameters for structural analysis (diagnostic): */ - double dvx[2]; - double dvy[2]; - double vx,vy; - double dvxdx,dvydy; - double v_gauss[2]={0.0}; - - double K[2][2]={0.0}; - double KDL[2][2]={0.0}; - int found=0; - int dim; - - /*inputs: */ - Input* vxaverage_input=NULL; - Input* vyaverage_input=NULL; - bool artdiff; + + /*Intermediaries */ + bool artdiff; + int i,j,ig,dim; + double Jdettria ,vx,vy,dvxdx,dvydy; + double dvx[2],dvy[2]; + double xyz_list[NUMVERTICES][3]; + double L[NUMVERTICES]; + double B[2][NUMVERTICES]; + double Bprime[2][NUMVERTICES]; + double K[2][2] = {0.0}; + double KDL[2][2] = {0.0}; + double DL[2][2] = {0.0}; + double DLprime[2][2] = {0.0}; + double DL_scalar; + double Ke_gg[numdof][numdof] = {0.0}; + double Ke_gg_gaussian[numdof][numdof] = {0.0}; + double Ke_gg_thickness1[numdof][numdof] = {0.0}; + double Ke_gg_thickness2[numdof][numdof] = {0.0}; + int *doflist = NULL; + GaussTria *gauss = NULL; /* Get node coordinates and dof list: */ @@ -2442,9 +2422,9 @@ GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - /*retrieve some parameters: */ + /*Retrieve all Inputs and parameters: */ this->parameters->FindParam(&artdiff,ArtDiffEnum); this->parameters->FindParam(&dim,DimEnum); - - /*Retrieve all inputs we will be needed*/ + Input* vxaverage_input=NULL; + Input* vyaverage_input=NULL; if(dim==2){ vxaverage_input=inputs->GetInput(VxEnum); ISSMASSERT(vxaverage_input); @@ -2458,5 +2438,4 @@ //Create Artificial diffusivity once for all if requested if(artdiff){ - //Get the Jacobian determinant gauss=new GaussTria(); gauss->GaussCenter(); @@ -2464,10 +2443,8 @@ delete gauss; - //Build K matrix (artificial diffusivity matrix) - vxaverage_input->GetParameterAverage(&v_gauss[0]); - vyaverage_input->GetParameterAverage(&v_gauss[1]); - - K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]); - K[1][1]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[1]); + vxaverage_input->GetParameterAverage(&vx); + vyaverage_input->GetParameterAverage(&vy); + K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(vx); + K[1][1]=pow(Jdettria,(double).5)/2.0*fabs(vy); } @@ -2478,15 +2455,10 @@ gauss->GaussPoint(ig); - /* Get Jacobian determinant: */ GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); - - /*Get B and B prime matrix: */ GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss); GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss); - //Get vx, vy and their derivatives at gauss point vxaverage_input->GetParameterValue(&vx,gauss); vyaverage_input->GetParameterValue(&vy,gauss); - vxaverage_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],gauss); vyaverage_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],gauss); @@ -2494,8 +2466,6 @@ dvxdx=dvx[0]; dvydy=dvy[1]; - DL_scalar=gauss->weight*Jdettria; - //Create DL and DLprime matrix DL[0][0]=DL_scalar*dvxdx; DL[1][1]=DL_scalar*dvydy; @@ -2503,7 +2473,4 @@ DLprime[0][0]=DL_scalar*vx; DLprime[1][1]=DL_scalar*vy; - - //Do the triple product tL*D*L. - //Ke_gg_thickness=B'*DLprime*Bprime; TripleMultiply( &B[0][0],2,numdof,1, @@ -2517,11 +2484,8 @@ &Ke_gg_thickness2[0][0],0); - /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */ for( i=0; iNewElementMatrix(MacAyealApproximationEnum); - - /*Add Ke_gg values to Ke element stifness matrix: */ if(Ke_gg_viscous) Ke->AddValues(Ke_gg_viscous); if(Ke_gg_friction)Ke->AddValues(Ke_gg_friction); @@ -2819,7 +2776,5 @@ xfree((void**)&Ke_gg_viscous); xfree((void**)&Ke_gg_friction); - -} - +} /*}}}*/ /*FUNCTION Tria::CreateKMatrixDiagnosticMacAyealViscous{{{1*/ @@ -3547,5 +3502,4 @@ double K[2][2]={0.0}; double KDL[2][2]={0.0}; - int found; int dim; @@ -3708,5 +3662,4 @@ /*input parameters for structural analysis (diagnostic): */ double vx,vy; - int found; int dim; @@ -3838,5 +3791,4 @@ int i,j; - int found=0; /* node data: */ @@ -5098,5 +5050,5 @@ void Tria::CreatePVectorThermalShelf( Vec pg){ - int i,found; + int i; const int numdof=NUMVERTICES*NDOF1; @@ -5186,5 +5138,5 @@ void Tria::CreatePVectorThermalSheet( Vec pg){ - int i,found; + int i; const int numdof=NUMVERTICES*NDOF1; @@ -6016,5 +5968,5 @@ int* sinternaldoflist=NULL; int* sexternaldoflist=NULL; - bool symmetric=true; + bool square=true; /*retrieve some parameters: */ @@ -6037,7 +5989,7 @@ } - /*Use symmetric constructor for ElementMatrix: */ - if(!kff) Ke=new ElementMatrix(gsize,symmetric,gexternaldoflist); - else Ke=new ElementMatrix(gsize,symmetric,finternaldoflist,fexternaldoflist,fsize,sinternaldoflist,sexternaldoflist,ssize); + /*Use square constructor for ElementMatrix: */ + if(!kff) Ke=new ElementMatrix(gsize,square,gexternaldoflist); + else Ke=new ElementMatrix(gsize,square,finternaldoflist,fexternaldoflist,fsize,sinternaldoflist,sexternaldoflist,ssize); /*Free ressources and return:*/ @@ -6081,5 +6033,5 @@ } - /*Use symmetric constructor for ElementVector: */ + /*constructor for ElementVector: */ if(!kff)pe=new ElementVector(gsize,gexternaldoflist); else pe=new ElementVector(gsize,finternaldoflist,fexternaldoflist,fsize); Index: /issm/trunk/src/c/objects/Numerics/ElementMatrix.cpp =================================================================== --- /issm/trunk/src/c/objects/Numerics/ElementMatrix.cpp (revision 5820) +++ /issm/trunk/src/c/objects/Numerics/ElementMatrix.cpp (revision 5821) @@ -25,7 +25,9 @@ this->nrows=0; this->ncols=0; - this->symmetric=false; + this->values=NULL; + this->square=false; + this->kff=false; + this->row_fsize=0; - this->values=NULL; this->row_finternaldoflist=NULL; this->row_fexternaldoflist=NULL; @@ -33,4 +35,5 @@ this->row_sinternaldoflist=NULL; this->row_sexternaldoflist=NULL; + this->col_fsize=0; this->col_finternaldoflist=NULL; @@ -42,12 +45,12 @@ } /*}}}*/ -/*FUNCTION ElementMatrix::ElementMatrix(int gsize,bool symmetric,int* in_gexternaldoflist){{{1*/ -ElementMatrix::ElementMatrix(int gsize,bool symmetric,int* in_gexternaldoflist){ - - if(symmetric=false)ISSMERROR(" calling symmetric constructor with false symmetric flag!"); +/*FUNCTION ElementMatrix::ElementMatrix(int gsize,bool square,int* in_gexternaldoflist){{{1*/ +ElementMatrix::ElementMatrix(int gsize,bool square,int* in_gexternaldoflist){ + + if(square=false)ISSMERROR(" calling square constructor with false square flag!"); this->nrows=gsize; this->ncols=gsize; - this->symmetric=true; + this->square=true; this->kff=false; @@ -59,5 +62,7 @@ memcpy(this->gexternaldoflist,in_gexternaldoflist,this->nrows*sizeof(int)); } - else this->gexternaldoflist=NULL; + else{ + this->gexternaldoflist=NULL; + } /*don't do rows and cols, because kff is false:*/ @@ -78,12 +83,12 @@ } /*}}}*/ -/*FUNCTION ElementMatrix::ElementMatrix(int gsize,bool symmetric,int* finternaldoflist,int* fexternaldoflist,int fsize,int* sinternaldoflist,int* sexternaldoflist,int ssize){{{1*/ -ElementMatrix::ElementMatrix(int gsize,bool symmetric,int* finternaldoflist,int* fexternaldoflist,int fsize,int* sinternaldoflist,int* sexternaldoflist,int ssize){ - - if(symmetric=false)ISSMERROR(" calling symmetric constructor with false symmetric flag!"); +/*FUNCTION ElementMatrix::ElementMatrix(int gsize,bool square,int* finternaldoflist,int* fexternaldoflist,int fsize,int* sinternaldoflist,int* sexternaldoflist,int ssize){{{1*/ +ElementMatrix::ElementMatrix(int gsize,bool square,int* finternaldoflist,int* fexternaldoflist,int fsize,int* sinternaldoflist,int* sexternaldoflist,int ssize){ + + if(square=false)ISSMERROR(" calling square constructor with false square flag!"); this->nrows=gsize; this->ncols=gsize; - this->symmetric=true; + this->square=true; this->kff=true; @@ -149,9 +154,7 @@ void ElementMatrix::AddValues(double* Ke_gg){ - int i,j; - if(Ke_gg){ - for (i=0;inrows;i++){ - for(j=0;jncols;j++){ + for (int i=0;inrows;i++){ + for(int j=0;jncols;j++){ *(this->values+this->ncols*i+j)+=*(Ke_gg+this->ncols*i+j); } @@ -166,5 +169,5 @@ double* internalvalues=NULL; - if(this->symmetric){ + if(this->square){ /*only use row dofs to add values into global matrices: */ @@ -207,5 +210,5 @@ } else{ - ISSMERROR(" unsymmetric AddToGlobal routine not support yet!"); + ISSMERROR(" non square matrix AddToGlobal routine not support yet!"); } @@ -219,5 +222,5 @@ printf(" nrows: %i\n",nrows); printf(" ncols: %i\n",ncols); - printf(" symmetric: %s\n",symmetric?"true":"false"); + printf(" square: %s\n",square?"true":"false"); printf(" kff: %s\n",kff?"true":"false"); @@ -225,7 +228,5 @@ for(i=0;i