Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5838) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5839) @@ -705,5 +705,7 @@ /*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){ - CreateKMatrixDiagnosticMacAyeal( Kgg,Kff,Kfs); + ElementMatrix* Ke=CreateKMatrixDiagnosticMacAyeal(); + Ke->AddToGlobal(Kgg,Kff,Kfs); + delete Ke; } else if (analysis_type==DiagnosticHutterAnalysisEnum){ @@ -714,18 +716,12 @@ } else if (analysis_type==PrognosticAnalysisEnum){ - if (GetElementType()==P1Enum) - CreateKMatrixPrognostic_CG( Kgg); - else if (GetElementType()==P1DGEnum) - CreateKMatrixPrognostic_DG( Kgg); - else - ISSMERROR("Element type %s not supported yet",EnumToString(GetElementType())); + ElementMatrix* Ke=CreateKMatrixPrognostic(); + Ke->AddToGlobal(Kgg,Kff,Kfs); + delete Ke; } else if (analysis_type==BalancedthicknessAnalysisEnum || analysis_type==AdjointBalancedthicknessAnalysisEnum){ - if (GetElementType()==P1Enum) - CreateKMatrixBalancedthickness_CG( Kgg); - else if (GetElementType()==P1DGEnum) - CreateKMatrixBalancedthickness_DG( Kgg); - else - ISSMERROR("Element type %s not supported yet",EnumToString(GetElementType())); + ElementMatrix* Ke=CreateKMatrixBalancedthickness(); + Ke->AddToGlobal(Kgg,Kff,Kfs); + delete Ke; } else if (analysis_type==BalancedvelocitiesAnalysisEnum){ @@ -2391,6 +2387,20 @@ /*Tria specific routines: */ +/*FUNCTION Tria::CreateKMatrixBalancedthickness {{{1*/ +ElementMatrix* Tria::CreateKMatrixBalancedthickness(void){ + + switch(GetElementType()){ + case P1Enum: + return CreateKMatrixBalancedthickness_CG(); + case P1DGEnum: + return CreateKMatrixBalancedthickness_DG(); + default: + ISSMERROR("Element type %s not supported yet",EnumToString(GetElementType())); + } + +} +/*}}}*/ /*FUNCTION Tria::CreateKMatrixBalancedthickness_CG {{{1*/ -void Tria::CreateKMatrixBalancedthickness_CG(Mat Kgg){ +ElementMatrix* Tria::CreateKMatrixBalancedthickness_CG(void){ /*Constants*/ @@ -2415,12 +2425,12 @@ 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: */ + /*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); - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - - /*Retrieve all Inputs and parameters: */ this->parameters->FindParam(&artdiff,ArtDiffEnum); this->parameters->FindParam(&dim,DimEnum); @@ -2484,6 +2494,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]; if(artdiff){ @@ -2496,31 +2506,22 @@ &Ke_gg_gaussian[0][0],0); - for( i=0; ivalues[i*numdof+j]+=Ke_gg_gaussian[i][j]; } } - MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES); - - /*Clean up*/ + /*Clean up and return*/ delete gauss; - xfree((void**)&doflist); + return Ke; } /*}}}*/ /*FUNCTION Tria::CreateKMatrixBalancedthickness_DG {{{1*/ -void Tria::CreateKMatrixBalancedthickness_DG(Mat Kgg){ +ElementMatrix* Tria::CreateKMatrixBalancedthickness_DG(void){ /* local declarations */ + const int numdof=NDOF1*NUMVERTICES; int i,j; - - /* node data: */ - const int numdof=NDOF1*NUMVERTICES; + int ig; double xyz_list[NUMVERTICES][3]; - int* doflist=NULL; - - /* gaussian points: */ - int ig; GaussTria *gauss=NULL; - - /* matrices: */ double B[2][NUMVERTICES]; double Bprime[2][NUMVERTICES]; @@ -2531,15 +2532,14 @@ double Ke_gg2[numdof][numdof]={0.0}; double Jdettria; - - /*input parameters for structural analysis (diagnostic): */ double vx,vy; int dim; - /* Get node coordinates and dof list: */ + /*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); - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); this->parameters->FindParam(&dim,DimEnum); - - /*Retrieve all inputs we will be needed*/ Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); @@ -2551,16 +2551,12 @@ gauss->GaussPoint(ig); - /* Get Jacobian determinant: */ - GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); - - /*Get B and B prime matrix: */ /*WARNING: B and Bprime are inverted compared to usual prognostic!!!!*/ GetBPrognostic(&Bprime[0][0], &xyz_list[0][0], gauss); GetBprimePrognostic(&B[0][0], &xyz_list[0][0], gauss); - //Get vx, vy and their derivatives at gauss point vx_input->GetParameterValue(&vx,gauss); vy_input->GetParameterValue(&vy,gauss); + GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); DL_scalar=-gauss->weight*Jdettria; @@ -2568,5 +2564,4 @@ DLprime[1][1]=DL_scalar*vy; - //Do the triple product tL*D*L. TripleMultiply( &B[0][0],2,numdof,1, &DLprime[0][0],2,2,0, @@ -2574,15 +2569,10 @@ &Ke_gg2[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_gg2[i][j]; + } /*Clean up and return*/ delete gauss; - xfree((void**)&doflist); + return Ke; } /*}}}*/ @@ -2753,5 +2743,5 @@ /*}}}*/ /*FUNCTION Tria::CreateKMatrixDiagnosticMacAyeal {{{1*/ -void Tria::CreateKMatrixDiagnosticMacAyeal(Mat Kgg,Mat Kff, Mat Kfs){ +ElementMatrix* Tria::CreateKMatrixDiagnosticMacAyeal(void){ /*compute all stiffness matrices for this element*/ @@ -2760,11 +2750,9 @@ ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); - /*Add Ke element stiffness matrix to global stiffness matrix: */ - Ke->AddToGlobal(Kgg,Kff,Kfs); - - /*Free ressources:*/ + /*clean-up and return*/ delete Ke1; delete Ke2; - delete Ke; + return Ke; + } /*}}}*/ @@ -2828,5 +2816,5 @@ &Ke_g[0][0],0); - for(i=0;ivalues[numdof*j+i]+=Ke_g[i][j]; + for(i=0;ivalues[i*numdof+j]+=Ke_g[i][j]; } @@ -2909,5 +2897,5 @@ &Ke_g[0][0],0); - for(i=0;ivalues[numdof*j+i]+=Ke_g[i][j]; + for(i=0;ivalues[i*numdof+j]+=Ke_g[i][j]; } @@ -3348,21 +3336,28 @@ } /*}}}*/ +/*FUNCTION Tria::CreateKMatrixPrognostic {{{1*/ +ElementMatrix* Tria::CreateKMatrixPrognostic(void){ + + switch(GetElementType()){ + case P1Enum: + return CreateKMatrixPrognostic_CG(); + case P1DGEnum: + return CreateKMatrixPrognostic_DG(); + default: + ISSMERROR("Element type %s not supported yet",EnumToString(GetElementType())); + } + +} +/*}}}*/ /*FUNCTION Tria::CreateKMatrixPrognostic_CG {{{1*/ -void Tria::CreateKMatrixPrognostic_CG(Mat Kgg){ +ElementMatrix* Tria::CreateKMatrixPrognostic_CG(void){ /* local declarations */ + const int numdof=NDOF1*NUMVERTICES; int i,j; - - /* node data: */ - const int numdof=NDOF1*NUMVERTICES; + int ig; double xyz_list[NUMVERTICES][3]; - int* doflist=NULL; int numberofdofspernode=1; - - /* gaussian points: */ - int ig; GaussTria *gauss=NULL; - - /* matrices: */ double L[NUMVERTICES]; double B[2][NUMVERTICES]; @@ -3376,6 +3371,4 @@ double Ke_gg_thickness2[numdof][numdof]={0.0}; double Jdettria; - - /*input parameters for structural analysis (diagnostic): */ double dvx[2]; double dvy[2]; @@ -3386,21 +3379,18 @@ double KDL[2][2]={0.0}; int dim; - - /*inputs: */ - Input* vxaverage_input=NULL; - Input* vyaverage_input=NULL; double dt; bool artdiff; - /*retrieve some parameters: */ + /*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(&dt,DtEnum); this->parameters->FindParam(&dim,DimEnum); this->parameters->FindParam(&artdiff,ArtDiffEnum); - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - - /*Retrieve all inputs we will be needing: */ + Input* vxaverage_input=NULL; + Input* vyaverage_input=NULL; if(dim==2){ vxaverage_input=inputs->GetInput(VxEnum); ISSMASSERT(vxaverage_input); @@ -3434,13 +3424,9 @@ gauss->GaussPoint(ig); - /* Get Jacobian determinant: */ + GetL(&L[0], &xyz_list[0][0], gauss,numberofdofspernode); GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); - /*Get L matrix: */ - GetL(&L[0], &xyz_list[0][0], gauss,numberofdofspernode); - DL_scalar=gauss->weight*Jdettria; - /* Do the triple product tL*D*L: */ TripleMultiply( &L[0],1,numdof,1, &DL_scalar,1,1,0, @@ -3448,29 +3434,19 @@ &Ke_gg_gaussian[0][0],0); - /*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); - dvxdx=dvx[0]; dvydy=dvy[1]; - DL_scalar=dt*gauss->weight*Jdettria; - //Create DL and DLprime matrix DL[0][0]=DL_scalar*dvxdx; DL[1][1]=DL_scalar*dvydy; - DLprime[0][0]=DL_scalar*vx; DLprime[1][1]=DL_scalar*vy; - - //Do the triple product tL*D*L. - //Ke_gg_thickness=B'*DL*B+B'*DLprime*Bprime; TripleMultiply( &B[0][0],2,numdof,1, @@ -3484,12 +3460,9 @@ &Ke_gg_thickness2[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_gaussian[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){ - - /* Compute artificial diffusivity */ KDL[0][0]=DL_scalar*K[0][0]; KDL[1][1]=DL_scalar*K[1][1]; @@ -3500,36 +3473,24 @@ &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; } /*}}}*/ /*FUNCTION Tria::CreateKMatrixPrognostic_DG {{{1*/ -void Tria::CreateKMatrixPrognostic_DG(Mat Kgg){ +ElementMatrix* Tria::CreateKMatrixPrognostic_DG(void){ /* local declarations */ int i,j; - - /* node data: */ const int numdof=NDOF1*NUMVERTICES; double xyz_list[NUMVERTICES][3]; int* doflist=NULL; int numberofdofspernode=1; - - /* gaussian points: */ int ig; GaussTria *gauss=NULL; - - /* matrices: */ double L[NUMVERTICES]; double B[2][NUMVERTICES]; @@ -3542,23 +3503,18 @@ double Ke_gg2[numdof][numdof]={0.0}; double Jdettria; - - /*input parameters for structural analysis (diagnostic): */ double vx,vy; int dim; - - /*inputs: */ + double dt; + + /*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(&dt,DtEnum); + this->parameters->FindParam(&dim,DimEnum); Input* vxaverage_input=NULL; Input* vyaverage_input=NULL; - double dt; - - /*retrieve some parameters: */ - this->parameters->FindParam(&dt,DtEnum); - 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*/ if(dim==2){ vxaverage_input=inputs->GetInput(VxEnum); ISSMASSERT(vxaverage_input); @@ -3576,13 +3532,10 @@ gauss->GaussPoint(ig); - /* Get Jacobian determinant: */ + + GetL(&L[0], &xyz_list[0][0], gauss,numberofdofspernode); GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); - /*Get L matrix: */ - GetL(&L[0], &xyz_list[0][0], gauss,numberofdofspernode); - DL_scalar=gauss->weight*Jdettria; - /* Do the triple product tL*D*L: */ TripleMultiply( &L[0],1,numdof,1, &DL_scalar,1,1,0, @@ -3595,8 +3548,6 @@ GetBprimePrognostic(&B[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); - DL_scalar=-dt*gauss->weight*Jdettria; @@ -3604,5 +3555,4 @@ DLprime[1][1]=DL_scalar*vy; - //Do the triple product tL*D*L. TripleMultiply( &B[0][0],2,numdof,1, &DLprime[0][0],2,2,0, @@ -3610,16 +3560,11 @@ &Ke_gg2[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_gg1[i][j]; + for(i=0;ivalues[i*numdof+j]+=Ke_gg2[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 5838) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 5839) @@ -121,8 +121,9 @@ /*}}}*/ /*Tria specific routines:{{{1*/ - void CreateKMatrixBalancedthickness_DG(Mat Kgg); - void CreateKMatrixBalancedthickness_CG(Mat Kgg); + ElementMatrix* CreateKMatrixBalancedthickness(void); + ElementMatrix* CreateKMatrixBalancedthickness_DG(void); + ElementMatrix* CreateKMatrixBalancedthickness_CG(void); void CreateKMatrixBalancedvelocities(Mat Kgg); - void CreateKMatrixDiagnosticMacAyeal(Mat Kgg,Mat Kff, Mat Kfs); + ElementMatrix* CreateKMatrixDiagnosticMacAyeal(void); ElementMatrix* CreateKMatrixDiagnosticMacAyealViscous(void); ElementMatrix* CreateKMatrixDiagnosticMacAyealFriction(void); @@ -133,6 +134,7 @@ void CreateKMatrixDiagnosticSurfaceVert(Mat Kgg); void CreateKMatrixMelting(Mat Kgg); - void CreateKMatrixPrognostic_CG(Mat Kgg); - void CreateKMatrixPrognostic_DG(Mat Kgg); + ElementMatrix* CreateKMatrixPrognostic(void); + ElementMatrix* CreateKMatrixPrognostic_CG(void); + ElementMatrix* CreateKMatrixPrognostic_DG(void); void CreateKMatrixSlope(Mat Kgg); void CreateKMatrixThermal(Mat Kgg); Index: /issm/trunk/src/c/objects/Numerics/ElementMatrix.cpp =================================================================== --- /issm/trunk/src/c/objects/Numerics/ElementMatrix.cpp (revision 5838) +++ /issm/trunk/src/c/objects/Numerics/ElementMatrix.cpp (revision 5839) @@ -358,14 +358,14 @@ for(i=0;i