Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5839) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5840) @@ -695,5 +695,6 @@ void Tria::CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs){ - /*retrive parameters: */ + /*retreive parameters: */ + ElementMatrix* Ke=NULL; int analysis_type; parameters->FindParam(&analysis_type,AnalysisTypeEnum); @@ -704,32 +705,32 @@ /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ - if (analysis_type==DiagnosticHorizAnalysisEnum || analysis_type==AdjointHorizAnalysisEnum){ - ElementMatrix* Ke=CreateKMatrixDiagnosticMacAyeal(); + switch(analysis_type){ + case DiagnosticHorizAnalysisEnum: case AdjointHorizAnalysisEnum: + Ke=CreateKMatrixDiagnosticMacAyeal(); + break; + case DiagnosticHutterAnalysisEnum: + Ke=CreateKMatrixDiagnosticHutter(); + break; + case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum: + Ke=CreateKMatrixSlope(); + break; + case PrognosticAnalysisEnum: + Ke=CreateKMatrixPrognostic(); + break; + case BalancedthicknessAnalysisEnum: case AdjointBalancedthicknessAnalysisEnum: + Ke=CreateKMatrixBalancedthickness(); + break; + case BalancedvelocitiesAnalysisEnum: + Ke=CreateKMatrixBalancedvelocities(); + break; + default: + ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type)); + } + + /*Add to global matrix*/ + if(Ke){ Ke->AddToGlobal(Kgg,Kff,Kfs); delete Ke; } - else if (analysis_type==DiagnosticHutterAnalysisEnum){ - CreateKMatrixDiagnosticHutter( Kgg); - } - else if (analysis_type==BedSlopeXAnalysisEnum || analysis_type==SurfaceSlopeXAnalysisEnum || analysis_type==BedSlopeYAnalysisEnum || analysis_type==SurfaceSlopeYAnalysisEnum){ - CreateKMatrixSlope( Kgg); - } - else if (analysis_type==PrognosticAnalysisEnum){ - ElementMatrix* Ke=CreateKMatrixPrognostic(); - Ke->AddToGlobal(Kgg,Kff,Kfs); - delete Ke; - } - else if (analysis_type==BalancedthicknessAnalysisEnum || analysis_type==AdjointBalancedthicknessAnalysisEnum){ - ElementMatrix* Ke=CreateKMatrixBalancedthickness(); - Ke->AddToGlobal(Kgg,Kff,Kfs); - delete Ke; - } - else if (analysis_type==BalancedvelocitiesAnalysisEnum){ - CreateKMatrixBalancedvelocities( Kgg); - } - else{ - ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type)); - } - } /*}}}*/ @@ -2494,6 +2495,6 @@ &Ke_gg_thickness2[0][0],0); - for( i=0; ivalues[i*numdof+j]+=Ke_gg_thickness1[i][j]; - for( i=0; ivalues[i*numdof+j]+=Ke_gg_thickness2[i][j]; + for(i=0;ivalues[i*numdof+j]+=Ke_gg_thickness1[i][j]; + for(i=0;ivalues[i*numdof+j]+=Ke_gg_thickness2[i][j]; if(artdiff){ @@ -2506,5 +2507,5 @@ &Ke_gg_gaussian[0][0],0); - for( i=0; ivalues[i*numdof+j]+=Ke_gg_gaussian[i][j]; + for(i=0;ivalues[i*numdof+j]+=Ke_gg_gaussian[i][j]; } } @@ -2578,19 +2579,12 @@ /*}}}*/ /*FUNCTION Tria::CreateKMatrixBalancedvelocities {{{1*/ -void Tria::CreateKMatrixBalancedvelocities(Mat Kgg){ +ElementMatrix* Tria::CreateKMatrixBalancedvelocities(void){ /* local declarations */ + const int numdof=NDOF1*NUMVERTICES; int i,j; - - /* node data: */ - 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]; @@ -2604,6 +2598,4 @@ double Ke_gg_velocities2[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point. double Jdettria; - - /*input parameters for structural analysis (diagnostic): */ double surface_normal[3]; double surface_list[3]; @@ -2617,20 +2609,17 @@ double KDL[2][2]={0.0}; int dim; - - /*inputs: */ + bool artdiff; + + /*Initialize Element matrix and return if necessary*/ + if(IsOnWater()) return NULL; + ElementMatrix* Ke=this->NewElementMatrix(NoneApproximationEnum); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + this->parameters->FindParam(&artdiff,ArtDiffEnum); + this->parameters->FindParam(&dim,DimEnum); + Input* surface_input=inputs->GetInput(SurfaceEnum); ISSMASSERT(surface_input); Input* vxaverage_input=NULL; Input* vyaverage_input=NULL; - bool artdiff; - - /*retrieve some parameters: */ - this->parameters->FindParam(&artdiff,ArtDiffEnum); - this->parameters->FindParam(&dim,DimEnum); - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - - /*Retrieve all inputs we will be needed*/ - Input* surface_input=inputs->GetInput(SurfaceEnum); ISSMASSERT(surface_input); if(dim==2){ vxaverage_input=inputs->GetInput(VxEnum); ISSMASSERT(vxaverage_input); @@ -2664,5 +2653,4 @@ //Create Artificial diffusivity once for all if requested if(artdiff){ - //Get the Jacobian determinant gauss=new GaussTria(); gauss->GaussCenter(); @@ -2670,5 +2658,4 @@ delete gauss; - //Build K matrix (artificial diffusivity matrix) vxaverage_input->GetParameterAverage(&v_gauss[0]); vyaverage_input->GetParameterAverage(&v_gauss[1]); @@ -2684,18 +2671,12 @@ 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 + GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); + 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); - dvxdx=dvx[0]; dvydy=dvy[1]; @@ -2706,6 +2687,4 @@ DLprime[1][1]=DL_scalar*ny; - //Do the triple product tL*D*L. - //Ke_gg_velocities=B'*DLprime*Bprime; TripleMultiply( &B[0][0],2,numdof,1, &DLprime[0][0],2,2,0, @@ -2713,10 +2692,7 @@ &Ke_gg_velocities2[0][0],0); - /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */ - for( i=0; ivalues[i*numdof+j]+=Ke_gg_velocities2[i][j]; if(artdiff){ - - /* Compute artificial diffusivity */ KDL[0][0]=DL_scalar*K[0][0]; KDL[1][1]=DL_scalar*K[1][1]; @@ -2727,17 +2703,13 @@ &Ke_gg_gaussian[0][0],0); - /* Add artificial diffusivity matrix */ - for( i=0; ivalues[i*numdof+j]+=Ke_gg_gaussian[i][j]; } } - - /*Add Ke_gg to global matrix Kgg: */ - MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES); /*Clean up and return*/ delete gauss; - xfree((void**)&doflist); + return Ke; } /*}}}*/ @@ -3139,29 +3111,23 @@ /*}}}*/ /*FUNCTION Tria::CreateKMatrixDiagnosticHutter{{{1*/ -void Tria::CreateKMatrixDiagnosticHutter(Mat Kgg){ - - /*Collapsed formulation: */ - int i; - int connectivity; - const int numdofs=NUMVERTICES*NDOF2; - int* doflist=NULL; - double Ke_gg[numdofs][numdofs]={0.0}; - - /*If on water, skip: */ - if(IsOnWater())return; - - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - - /*Spawn 3 sing elements: */ - for(i=0;i<3;i++){ +ElementMatrix* Tria::CreateKMatrixDiagnosticHutter(void){ + + /*Intermediaries*/ + const int numdof=NUMVERTICES*NDOF2; + int i,connectivity; + + /*Initialize Element matrix and return if necessary*/ + if(IsOnWater()) return NULL; + ElementMatrix* Ke=this->NewElementMatrix(NoneApproximationEnum); + + /*Create Element matrix*/ + for(i=0;iGetConnectivity(); - Ke_gg[2*i][2*i]=1/(double)connectivity; - Ke_gg[2*i+1][2*i+1]=1/(double)connectivity; - } - - MatSetValues(Kgg,numdofs,doflist,numdofs,doflist,(const double*)Ke_gg,ADD_VALUES); - - /*Free ressources:*/ - xfree((void**)&doflist); + Ke->values[(2*i)*numdof +(2*i) ]=1/(double)connectivity; + Ke->values[(2*i+1)*numdof+(2*i+1)]=1/(double)connectivity; + } + + /*Clean up and return*/ + return Ke; } /*}}}*/ @@ -3570,22 +3536,22 @@ /*}}}*/ /*FUNCTION Tria::CreateKMatrixSlope {{{1*/ -void Tria::CreateKMatrixSlope(Mat Kgg){ +ElementMatrix* Tria::CreateKMatrixSlope(void){ /*constants: */ - const int numnodes=3; - const int numdof=NDOF1*numnodes; + const int numdof=NDOF1*NUMVERTICES; /* Intermediaries */ int i,j,ig; double DL_scalar,Jdet; - double xyz_list[numnodes][3]; + double xyz_list[NUMVERTICES][3]; double L[1][3]; - double Ke[numdof][numdof] = {0.0}; double Ke_g[numdof][numdof]; GaussTria *gauss = NULL; - int *doflist = NULL; - - GetVerticesCoordinates(&xyz_list[0][0], nodes, numnodes); - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + + /*Initialize Element matrix and return if necessary*/ + if(IsOnWater()) return NULL; + ElementMatrix* Ke=this->NewElementMatrix(NoneApproximationEnum); + + GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES); /* Start looping on the number of gaussian points: */ @@ -3605,12 +3571,10 @@ &Ke_g[0][0],0); - for(i=0;ivalues[i*numdof+j]+=Ke_g[i][j]; + } + + /*Clean up and return*/ delete gauss; - xfree((void**)&doflist); + return Ke; } /*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 5839) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 5840) @@ -124,5 +124,5 @@ ElementMatrix* CreateKMatrixBalancedthickness_DG(void); ElementMatrix* CreateKMatrixBalancedthickness_CG(void); - void CreateKMatrixBalancedvelocities(Mat Kgg); + ElementMatrix* CreateKMatrixBalancedvelocities(void); ElementMatrix* CreateKMatrixDiagnosticMacAyeal(void); ElementMatrix* CreateKMatrixDiagnosticMacAyealViscous(void); @@ -131,5 +131,5 @@ void CreateKMatrixCouplingMacAyealPattynFriction(Mat Kgg); void CreateKMatrixDiagnosticPattynFriction(Mat Kgg); - void CreateKMatrixDiagnosticHutter(Mat Kgg); + ElementMatrix* CreateKMatrixDiagnosticHutter(void); void CreateKMatrixDiagnosticSurfaceVert(Mat Kgg); void CreateKMatrixMelting(Mat Kgg); @@ -137,5 +137,5 @@ ElementMatrix* CreateKMatrixPrognostic_CG(void); ElementMatrix* CreateKMatrixPrognostic_DG(void); - void CreateKMatrixSlope(Mat Kgg); + ElementMatrix* CreateKMatrixSlope(void); void CreateKMatrixThermal(Mat Kgg); void CreatePVectorBalancedthickness_DG(Vec pg);