Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5830) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5831) @@ -2755,17 +2755,8 @@ void Tria::CreateKMatrixDiagnosticMacAyeal(Mat Kgg,Mat Kff, Mat Kfs){ - /*Intermediaries*/ - double *Ke_gg_viscous = NULL; - double *Ke_gg_friction = NULL; - ElementMatrix *Ke = NULL; - - /*compute stiffness matrices on the g-set, at the element level: */ - Ke_gg_viscous = CreateKMatrixDiagnosticMacAyealViscous(); - Ke_gg_friction = CreateKMatrixDiagnosticMacAyealFriction(); - - /*Add Ke_gg values to Ke element stifness matrix: */ - Ke=this->NewElementMatrix(MacAyealApproximationEnum); - if(Ke_gg_viscous) Ke->AddValues(Ke_gg_viscous); - if(Ke_gg_friction)Ke->AddValues(Ke_gg_friction); + /*compute all stiffness matrices for this element*/ + ElementMatrix* Ke1=CreateKMatrixDiagnosticMacAyealViscous(); + ElementMatrix* Ke2=CreateKMatrixDiagnosticMacAyealFriction(); + ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); /*Add Ke element stiffness matrix to global stiffness matrix: */ @@ -2773,66 +2764,36 @@ /*Free ressources:*/ + delete Ke1; + delete Ke2; delete Ke; - xfree((void**)&Ke_gg_viscous); - xfree((void**)&Ke_gg_friction); } /*}}}*/ /*FUNCTION Tria::CreateKMatrixDiagnosticMacAyealViscous{{{1*/ -double* Tria::CreateKMatrixDiagnosticMacAyealViscous(void){ - - /* local declarations */ - int i,j; - - /*output: */ - double* Ke_gg=NULL; - - /* node data: */ - const int numdof=2*NUMVERTICES; - double xyz_list[NUMVERTICES][3]; - - /* gaussian points: */ - int ig; - GaussTria *gauss=NULL; - - /* material data: */ - double viscosity; //viscosity - double newviscosity; //viscosity - double oldviscosity; //viscosity - - /* strain rate: */ - double epsilon[3]; /* epsilon=[exx,eyy,exy];*/ - double oldepsilon[3]; /* oldepsilon=[exx,eyy,exy];*/ - - /* matrices: */ - double B[3][numdof]; - double Bprime[3][numdof]; - double D[3][3]={0.0}; // material matrix, simple scalar matrix. - double D_scalar; - - /*parameters: */ - double viscosity_overshoot; - - /* local element matrices: */ - double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point. - - double Jdet; - - /*input parameters for structural analysis (diagnostic): */ - double thickness; - - - /*retrieve some parameters: */ - this->parameters->FindParam(&viscosity_overshoot,ViscosityOvershootEnum); - - /*First, if we are on water, return empty matrix: */ - if(IsOnWater()) return Ke_gg; - - /*allocate output: */ - Ke_gg=(double*)xcalloc(numdof*numdof,sizeof(double)); - - /* Get node coordinates and dof list: */ +ElementMatrix* Tria::CreateKMatrixDiagnosticMacAyealViscous(void){ + + /*Constants*/ + const int numdof=NDOF2*NUMVERTICES; + + /*Intermediaries*/ + int i,j,ig; + double xyz_list[NUMVERTICES][3]; + GaussTria *gauss = NULL; + double viscosity,newviscosity,oldviscosity; + double viscosity_overshoot,thickness; + double epsilon[3]; /* epsilon=[exx,eyy,exy]; */ + double oldepsilon[3]; /* oldepsilon=[exx,eyy,exy]; */ + double B[3][numdof]; + double Bprime[3][numdof]; + double D[3][3] = {0.0}; + double D_scalar; + double Ke_g[numdof][numdof]; + double Jdet; + + /*Initialize Element matrix and return if necessary*/ + if(IsOnWater()) return NULL; + ElementMatrix* Ke=this->NewElementMatrix(MacAyealApproximationEnum); + + /*Retrieve all inputs and parameters*/ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - - /*Retrieve all inputs we will be needing: */ Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); @@ -2840,4 +2801,5 @@ Input* vxold_input=inputs->GetInput(VxOldEnum); ISSMASSERT(vxold_input); Input* vyold_input=inputs->GetInput(VyOldEnum); ISSMASSERT(vyold_input); + this->parameters->FindParam(&viscosity_overshoot,ViscosityOvershootEnum); /* Start looping on the number of gaussian points: */ @@ -2847,135 +2809,79 @@ gauss->GaussPoint(ig); - /*Compute thickness at gaussian point: */ thickness_input->GetParameterValue(&thickness, gauss); - - /*Get strain rate from velocity: */ this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); this->GetStrainRate2d(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input); - - /*Get viscosity: */ matice->GetViscosity2d(&viscosity, &epsilon[0]); matice->GetViscosity2d(&oldviscosity, &oldepsilon[0]); - - /* Get Jacobian determinant: */ GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); - /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant - onto this scalar matrix, so that we win some computational time: */ newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); D_scalar=2*newviscosity*thickness*gauss->weight*Jdet; - - for (i=0;i<3;i++){ - D[i][i]=D_scalar; - } - - /*Get B and Bprime matrices: */ + for (i=0;i<3;i++) D[i][i]=D_scalar; + GetBMacAyeal(&B[0][0], &xyz_list[0][0], gauss); GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss); - /* Do the triple product tB*D*Bprime: */ - TripleMultiply( &B[0][0],3,numdof,1, + TripleMultiply(&B[0][0],3,numdof,1, &D[0][0],3,3,0, &Bprime[0][0],3,numdof,0, - &Ke_gg_gaussian[0][0],0); - - /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */ - for( i=0; ivalues[numdof*j+i]+=Ke_g[i][j]; } /*Clean up and return*/ delete gauss; - return Ke_gg; + return Ke; } /*}}}*/ /*FUNCTION Tria::CreateKMatrixDiagnosticMacAyealFriction {{{1*/ void Tria::CreateKMatrixDiagnosticMacAyealFriction(Mat Kgg,Mat Kff, Mat Kfs){ - - double *Ke_gg_friction = NULL; - ElementMatrix *Ke = NULL; - - /*compute stiffness matrices on the g-set, at the element level: */ - Ke_gg_friction = CreateKMatrixDiagnosticMacAyealFriction(); - - if(Ke_gg_friction){ - Ke=this->NewElementMatrix(MacAyealApproximationEnum); - Ke->AddValues(Ke_gg_friction); + + ElementMatrix* Ke=CreateKMatrixDiagnosticMacAyealFriction(); + + if(Ke){ Ke->AddToGlobal(Kgg,Kff,Kfs); delete Ke; } - - /*Free ressources:*/ - xfree((void**)&Ke_gg_friction); -} - +} /*}}}*/ /*FUNCTION Tria::CreateKMatrixDiagnosticMacAyealFriction {{{1*/ -double* Tria::CreateKMatrixDiagnosticMacAyealFriction(void){ - - /* local declarations */ - int i,j; - int analysis_type; - - /*output: */ - double* Ke_gg=NULL; - - /* node data: */ - const int numdof = 2 *NUMVERTICES; - double xyz_list[NUMVERTICES][3]; - int numberofdofspernode=2; - - /* gaussian points: */ - int ig; - GaussTria *gauss=NULL; - - /* matrices: */ - double L[2][numdof]; - double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag - double DL_scalar; - - /* local element matrices: */ - double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag - - double Jdet; - - /*slope: */ - double slope[2]={0.0,0.0}; - double slope_magnitude; - - /*friction: */ - Friction *friction = NULL; - double alpha2; - - double MAXSLOPE=.06; // 6 % - double MOUNTAINKEXPONENT=10; - - /*inputs: */ - int drag_type; - - - /*retrive parameters: */ - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - - /*retrieve inputs :*/ - inputs->GetParameterValue(&drag_type,DragTypeEnum); +ElementMatrix* Tria::CreateKMatrixDiagnosticMacAyealFriction(void){ + + /*Constants*/ + const int numdof=NDOF2*NUMVERTICES; + + /*Intermediaries*/ + int i,j,ig; + int analysis_type,drag_type; + double MAXSLOPE = .06; // 6 % + double MOUNTAINKEXPONENT = 10; + double slope_magnitude,alpha2; + double Jdet; + double L[2][numdof]; + double DL[2][2] = {{ 0,0 },{0,0}}; + double Ke_g[numdof][numdof]; + double DL_scalar; + double slope[2] = {0.0,0.0}; + double xyz_list[NUMVERTICES][3]; + Friction *friction = NULL; + GaussTria *gauss = NULL; + + /*Initialize Element matrix and return if necessary*/ + if(IsOnWater() || IsOnShelf()) return NULL; + ElementMatrix* Ke=this->NewElementMatrix(MacAyealApproximationEnum); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); Input* surface_input=inputs->GetInput(SurfaceEnum); ISSMASSERT(surface_input); Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); Input* vz_input=inputs->GetInput(VzEnum); ISSMASSERT(vz_input); - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - - if (IsOnShelf()){ - /*no friction, do nothing*/ - return Ke_gg; - } - - /*allocte Ke_gg matrix: */ - Ke_gg=(double*)xcalloc(numdof*numdof,sizeof(double)); + inputs->GetParameterValue(&drag_type,DragTypeEnum); + parameters->FindParam(&analysis_type,AnalysisTypeEnum); /*build friction object, used later on: */ - if (drag_type!=2)ISSMERROR(" non-viscous friction not supported yet!"); + if (drag_type!=2) ISSMERROR(" non-viscous friction not supported yet!"); friction=new Friction("2d",inputs,matpar,analysis_type); @@ -2985,7 +2891,4 @@ gauss->GaussPoint(ig); - - /*Friction: */ - friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); // If we have a slope > 6% for this element, it means we are on a mountain. In this particular case, @@ -2993,29 +2896,18 @@ surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss); slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2)); - - if (slope_magnitude>MAXSLOPE){ - alpha2=pow((double)10,MOUNTAINKEXPONENT); - } - - /* Get Jacobian determinant: */ + if(slope_magnitude>MAXSLOPE) alpha2=pow((double)10,MOUNTAINKEXPONENT); + else friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); + + GetL(&L[0][0], &xyz_list[0][0], gauss,NDOF2); GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); - - /*Get L matrix: */ - GetL(&L[0][0], &xyz_list[0][0], gauss,numberofdofspernode); - + DL_scalar=alpha2*gauss->weight*Jdet; + for (i=0;i<2;i++) DL[i][i]=DL_scalar; - DL_scalar=alpha2*gauss->weight*Jdet; - for (i=0;i<2;i++){ - DL[i][i]=DL_scalar; - } - - /* Do the triple producte tL*D*L: */ TripleMultiply( &L[0][0],2,numdof,1, &DL[0][0],2,2,0, &L[0][0],2,numdof,0, - &Ke_gg_gaussian[0][0],0); - - for( i=0; ivalues[numdof*j+i]+=Ke_g[i][j]; } @@ -3023,5 +2915,5 @@ delete gauss; delete friction; - return Ke_gg; + return Ke; } /*}}}*/ @@ -5954,9 +5846,9 @@ int fsize; int ssize; - int* gexternaldoflist=NULL; - int* finternaldoflist=NULL; - int* fexternaldoflist=NULL; - int* sinternaldoflist=NULL; - int* sexternaldoflist=NULL; + int* gglobaldoflist=NULL; + int* flocaldoflist=NULL; + int* fglobaldoflist=NULL; + int* slocaldoflist=NULL; + int* sglobaldoflist=NULL; bool square=true; @@ -5972,22 +5864,22 @@ /*get dof lists for f and s set: */ - if(!kff) gexternaldoflist=this->GetExternalDofList(approximation,GsetEnum); - else{ - finternaldoflist=this->GetInternalDofList(approximation,FsetEnum); - fexternaldoflist=this->GetExternalDofList(approximation,FsetEnum); - sinternaldoflist=this->GetInternalDofList(approximation,SsetEnum); - sexternaldoflist=this->GetExternalDofList(approximation,SsetEnum); + gglobaldoflist=this->GetExternalDofList(approximation,GsetEnum); + if(kff){ + flocaldoflist=this->GetInternalDofList(approximation,FsetEnum); + fglobaldoflist=this->GetExternalDofList(approximation,FsetEnum); + slocaldoflist=this->GetInternalDofList(approximation,SsetEnum); + sglobaldoflist=this->GetExternalDofList(approximation,SsetEnum); } /*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); + if(!kff) Ke=new ElementMatrix(square,gglobaldoflist,gsize); + else Ke=new ElementMatrix(square,gglobaldoflist,gsize,flocaldoflist,fglobaldoflist,fsize,slocaldoflist,sglobaldoflist,ssize); /*Free ressources and return:*/ - xfree((void**)&gexternaldoflist); - xfree((void**)&finternaldoflist); - xfree((void**)&fexternaldoflist); - xfree((void**)&sinternaldoflist); - xfree((void**)&sexternaldoflist); + xfree((void**)&gglobaldoflist); + xfree((void**)&flocaldoflist); + xfree((void**)&fglobaldoflist); + xfree((void**)&slocaldoflist); + xfree((void**)&sglobaldoflist); return Ke; @@ -6004,7 +5896,7 @@ int gsize; int fsize; - int* gexternaldoflist=NULL; - int* finternaldoflist=NULL; - int* fexternaldoflist=NULL; + int* gglobaldoflist=NULL; + int* flocaldoflist=NULL; + int* fglobaldoflist=NULL; /*retrieve some parameters: */ @@ -6017,19 +5909,19 @@ /*get dof lists for f and s set: */ if(!kff){ - gexternaldoflist=this->GetExternalDofList(approximation,GsetEnum); + gglobaldoflist=this->GetExternalDofList(approximation,GsetEnum); } else{ - finternaldoflist=this->GetInternalDofList(approximation,FsetEnum); - fexternaldoflist=this->GetExternalDofList(approximation,FsetEnum); + flocaldoflist=this->GetInternalDofList(approximation,FsetEnum); + fglobaldoflist=this->GetExternalDofList(approximation,FsetEnum); } /*constructor for ElementVector: */ - if(!kff)pe=new ElementVector(gsize,gexternaldoflist); - else pe=new ElementVector(gsize,finternaldoflist,fexternaldoflist,fsize); + if(!kff)pe=new ElementVector(gsize,gglobaldoflist); + else pe=new ElementVector(gsize,flocaldoflist,fglobaldoflist,fsize); /*Free ressources and return:*/ - xfree((void**)&gexternaldoflist); - xfree((void**)&finternaldoflist); - xfree((void**)&fexternaldoflist); + xfree((void**)&gglobaldoflist); + xfree((void**)&flocaldoflist); + xfree((void**)&fglobaldoflist); return pe; Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 5830) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 5831) @@ -125,7 +125,7 @@ void CreateKMatrixBalancedvelocities(Mat Kgg); void CreateKMatrixDiagnosticMacAyeal(Mat Kgg,Mat Kff, Mat Kfs); - double* CreateKMatrixDiagnosticMacAyealViscous(void); - double* CreateKMatrixDiagnosticMacAyealFriction(void); - void CreateKMatrixDiagnosticMacAyealFriction(Mat Kgg,Mat Kff, Mat Kfs); + ElementMatrix* CreateKMatrixDiagnosticMacAyealViscous(void); + ElementMatrix* CreateKMatrixDiagnosticMacAyealFriction(void); + void CreateKMatrixDiagnosticMacAyealFriction(Mat Kgg,Mat Kff, Mat Kfs); void CreateKMatrixCouplingMacAyealPattynFriction(Mat Kgg); void CreateKMatrixDiagnosticPattynFriction(Mat Kgg); Index: /issm/trunk/src/c/objects/Numerics/ElementMatrix.cpp =================================================================== --- /issm/trunk/src/c/objects/Numerics/ElementMatrix.cpp (revision 5830) +++ /issm/trunk/src/c/objects/Numerics/ElementMatrix.cpp (revision 5831) @@ -45,6 +45,121 @@ } /*}}}*/ -/*FUNCTION ElementMatrix::ElementMatrix(int gsize,bool square,int* in_gglobaldoflist){{{1*/ -ElementMatrix::ElementMatrix(int gsize,bool square,int* in_gglobaldoflist){ +/*FUNCTION ElementMatrix::ElementMatrix(ElementMatrix* Ke1, ElementMatrix* Ke2){{{1*/ +ElementMatrix::ElementMatrix(ElementMatrix* Ke1, ElementMatrix* Ke2){ + + /*intermediaries*/ + int i,j,counter; + int gsize,fsize,ssize; + int* P=NULL; + bool found; + + /*If one of the two matrix is NULL, we copy the other one*/ + if(!Ke1 && !Ke2){ + ISSMERROR("Two input element matrices are NULL"); + } + else if(!Ke1){ + this->Init(Ke2); + return; + } + else if(!Ke2){ + this->Init(Ke1); + return; + } + + /*General Case: Ke1 and Ke2 are not empty*/ + if(!Ke1->square || !Ke2->square) ISSMERROR("merging 2 non square matrices not implemented yet"); + + /*Initialize itransformation matrix Ke[P[i]] = Ke2[i]*/ + P=(int*)xmalloc(Ke2->nrows*sizeof(int)); + + /*1: Get the new numbering of Ke2 and get size of the new matrix*/ + gsize=Ke1->nrows; + for(i=0;inrows;i++){ + found=false; + for(j=0;jnrows;j++){ + if(Ke2->gglobaldoflist[i]==Ke1->gglobaldoflist[j]){ + found=true; P[i]=j; break; + } + } + if(!found){ + P[i]=gsize; gsize++; + } + } + + /*2: Initialize static fields*/ + this->nrows=gsize; + this->ncols=gsize; + this->square=true; + this->kff=Ke1->kff; + + /*Gset and values*/ + this->gglobaldoflist=(int*)xmalloc(this->nrows*sizeof(int)); + this->values=(double*)xcalloc(this->nrows*this->ncols,sizeof(double)); + for(i=0;inrows;i++){ + for(j=0;jncols;j++){ + this->values[i*this->ncols+j] += Ke1->values[i*Ke1->ncols+j]; + } + this->gglobaldoflist[i]=Ke1->gglobaldoflist[i]; + } + for(i=0;inrows;i++){ + for(j=0;jncols;j++){ + this->values[P[i]*this->ncols+P[j]] += Ke2->values[i*Ke2->ncols+j]; + } + this->gglobaldoflist[P[i]]=Ke2->gglobaldoflist[i]; + } + + /*Fset*/ + fsize=Ke1->row_fsize; + for(i=0;irow_fsize;i++){ + if(P[Ke2->row_flocaldoflist[i]] >= Ke1->nrows) fsize++; + } + printf("fsize = %i\n",fsize); + ISSMERROR("STOP"); + if(fsize){ + this->row_flocaldoflist =(int*)xmalloc(fsize*sizeof(int)); + this->row_fglobaldoflist=(int*)xmalloc(fsize*sizeof(int)); + for(i=0;irow_fsize;i++){ + this->row_flocaldoflist[i]=Ke1->row_flocaldoflist[i]; + } + counter=Ke1->row_fsize; + for(i=0;irow_fsize;i++){ + if(P[Ke2->row_flocaldoflist[i]] >= Ke1->nrows){ + this->row_flocaldoflist[counter]=P[Ke2->row_flocaldoflist[i]]; + } + } + } + else{ + this->row_flocaldoflist=NULL; + this->row_fglobaldoflist=NULL; + } + + ssize=Ke1->row_ssize+Ke2->row_ssize; + this->row_ssize=ssize; + if(ssize){ + ISSMERROR("STOP"); + this->row_slocaldoflist=(int*)xmalloc(ssize*sizeof(int)); + this->row_sglobaldoflist=(int*)xmalloc(ssize*sizeof(int)); + //memcpy(this->row_slocaldoflist,slocaldoflist,ssize*sizeof(int)); + //memcpy(this->row_sglobaldoflist,sglobaldoflist,ssize*sizeof(int)); + } + else{ + this->row_slocaldoflist=NULL; + this->row_sglobaldoflist=NULL; + } + + /*don't do cols, we can pick them up from the rows: */ + this->col_fsize=0; + this->col_flocaldoflist=NULL; + this->col_fglobaldoflist=NULL; + this->col_ssize=0; + this->col_slocaldoflist=NULL; + this->col_sglobaldoflist=NULL; + + /*clean-up*/ + xfree((void**)&P); +} +/*}}}*/ +/*FUNCTION ElementMatrix::ElementMatrix(bool square,int* in_gglobaldoflist,int gsize){{{1*/ +ElementMatrix::ElementMatrix(bool square,int* in_gglobaldoflist,int gsize){ if(square=false)ISSMERROR(" calling square constructor with false square flag!"); @@ -83,8 +198,9 @@ } /*}}}*/ -/*FUNCTION ElementMatrix::ElementMatrix(int gsize,bool square,int* flocaldoflist,int* fglobaldoflist,int fsize,int* slocaldoflist,int* sglobaldoflist,int ssize){{{1*/ -ElementMatrix::ElementMatrix(int gsize,bool square,int* flocaldoflist,int* fglobaldoflist,int fsize,int* slocaldoflist,int* sglobaldoflist,int ssize){ +/*FUNCTION ElementMatrix::ElementMatrix(bool square,int* in_gglobaldoflist,int gsize,int* flocaldoflist,int* fglobaldoflist,int fsize,int* slocaldoflist,int* sglobaldoflist,int ssize){{{1*/ +ElementMatrix::ElementMatrix(bool square,int* in_gglobaldoflist,int gsize,int* flocaldoflist,int* fglobaldoflist,int fsize,int* slocaldoflist,int* sglobaldoflist,int ssize){ if(square=false)ISSMERROR(" calling square constructor with false square flag!"); + ISSMASSERT(gsize); this->nrows=gsize; @@ -95,5 +211,9 @@ /*fill values with 0: */ this->values=(double*)xcalloc(this->nrows*this->ncols,sizeof(double)); - + + /*g dofs*/ + this->gglobaldoflist=(int*)xmalloc(gsize*sizeof(int)); + memcpy(this->gglobaldoflist,in_gglobaldoflist,gsize*sizeof(int)); + /*row dofs: */ this->row_fsize=fsize; @@ -121,5 +241,5 @@ } - /*don't do cols, we can't pick them up from the rows: */ + /*don't do cols, we can pick them up from the rows: */ this->col_fsize=0; this->col_flocaldoflist=NULL; @@ -128,7 +248,4 @@ this->col_slocaldoflist=NULL; this->col_sglobaldoflist=NULL; - - /*don't do g-set: */ - this->gglobaldoflist=NULL; } @@ -151,5 +268,5 @@ /*ElementMatrix specific routines: */ -/*FUNCTION ElementMatrix::AddValues(double* Ke_gg){{{1*/ +/*FUNCTION ElementMatrix::AddValues{{{1*/ void ElementMatrix::AddValues(double* Ke_gg){ @@ -163,5 +280,5 @@ } /*}}}*/ -/*FUNCTION ElementMatrix::AddToGlobal(Mat Kgg, Mat Kff, Mat Kfs){{{1*/ +/*FUNCTION ElementMatrix::AddToGlobal{{{1*/ void ElementMatrix::AddToGlobal(Mat Kgg, Mat Kff, Mat Kfs){ @@ -215,5 +332,5 @@ } /*}}}*/ -/*FUNCTION ElementMatrix::Echo(void){{{1*/ +/*FUNCTION ElementMatrix::Echo{{{1*/ void ElementMatrix::Echo(void){ @@ -262,2 +379,63 @@ } /*}}}*/ +/*FUNCTION ElementMatrix::Init{{{1*/ +void ElementMatrix::Init(ElementMatrix* Ke){ + + ISSMASSERT(Ke); + + this->nrows =Ke->nrows; + this->ncols =Ke->ncols; + this->square=Ke->square; + this->kff =Ke->kff; + this->values=(double*)xmalloc(this->nrows*this->ncols*sizeof(double)); + memcpy(this->values,Ke->values,this->nrows*this->ncols*sizeof(double)); + + this->row_fsize=Ke->row_fsize; + if(this->row_fsize){ + this->row_flocaldoflist=(int*)xmalloc(this->row_fsize*sizeof(int)); + memcpy(this->row_flocaldoflist,Ke->row_flocaldoflist,this->row_fsize*sizeof(int)); + this->row_fglobaldoflist=(int*)xmalloc(this->row_fsize*sizeof(int)); + memcpy(this->row_fglobaldoflist,Ke->row_fglobaldoflist,this->row_fsize*sizeof(int)); + } + else{ + this->row_flocaldoflist=NULL; + this->row_fglobaldoflist=NULL; + } + + this->row_ssize=Ke->row_ssize; + if(this->row_ssize){ + this->row_slocaldoflist=(int*)xmalloc(this->row_ssize*sizeof(int)); + memcpy(this->row_slocaldoflist,Ke->row_slocaldoflist,this->row_ssize*sizeof(int)); + this->row_sglobaldoflist=(int*)xmalloc(this->row_ssize*sizeof(int)); + memcpy(this->row_sglobaldoflist,Ke->row_sglobaldoflist,this->row_ssize*sizeof(int)); + } + else{ + this->row_slocaldoflist=NULL; + this->row_sglobaldoflist=NULL; + } + + this->col_fsize=Ke->col_fsize; + if(this->col_fsize){ + this->col_flocaldoflist=(int*)xmalloc(this->col_fsize*sizeof(int)); + memcpy(this->col_flocaldoflist,Ke->col_flocaldoflist,this->col_fsize*sizeof(int)); + this->col_fglobaldoflist=(int*)xmalloc(this->col_fsize*sizeof(int)); + memcpy(this->col_fglobaldoflist,Ke->col_fglobaldoflist,this->col_fsize*sizeof(int)); + } + else{ + this->col_flocaldoflist=NULL; + this->col_fglobaldoflist=NULL; + } + + this->col_ssize=Ke->col_ssize; + if(this->col_ssize){ + this->col_slocaldoflist=(int*)xmalloc(this->col_ssize*sizeof(int)); + memcpy(this->col_slocaldoflist,Ke->col_slocaldoflist,this->col_ssize*sizeof(int)); + this->col_sglobaldoflist=(int*)xmalloc(this->col_ssize*sizeof(int)); + memcpy(this->col_sglobaldoflist,Ke->col_sglobaldoflist,this->col_ssize*sizeof(int)); + } + else{ + this->col_slocaldoflist=NULL; + this->col_sglobaldoflist=NULL; + } +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Numerics/ElementMatrix.h =================================================================== --- /issm/trunk/src/c/objects/Numerics/ElementMatrix.h (revision 5830) +++ /issm/trunk/src/c/objects/Numerics/ElementMatrix.h (revision 5831) @@ -21,7 +21,7 @@ int nrows; int ncols; - double* values; bool square; bool kff; + double* values; //gset @@ -50,6 +50,7 @@ /*ElementMatrix constructors, destructors {{{1*/ ElementMatrix(); - ElementMatrix(int gsize,bool square,int* gglobaldoflist); - ElementMatrix(int gsize,bool square,int* flocaldoflist,int* fglobaldoflist,int fsize,int* slocaldoflist,int* sglobaldoflist,int ssize); + ElementMatrix(ElementMatrix* Ke1,ElementMatrix* Ke2); + ElementMatrix(bool square,int* gglobaldoflist,int gsize); + ElementMatrix(bool square,int* gglobaldoflist,int gsize,int* flocaldoflist,int* fglobaldoflist,int fsize,int* slocaldoflist,int* sglobaldoflist,int ssize); ~ElementMatrix(); /*}}}*/ @@ -58,4 +59,5 @@ void AddToGlobal(Mat Kgg, Mat Kff, Mat Kfs); void Echo(void); + void Init(ElementMatrix* Ke); /*}}}*/ };