Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5878) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5879) @@ -700,47 +700,37 @@ case DiagnosticHorizAnalysisEnum: case AdjointHorizAnalysisEnum: Ke=CreateKMatrixDiagnosticHoriz(); - if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL); - delete Ke; break; case DiagnosticHutterAnalysisEnum: Ke=CreateKMatrixDiagnosticHutter(); - if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL); - delete Ke; break; case DiagnosticVertAnalysisEnum: Ke=CreateKMatrixDiagnosticVert(); - if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL); - delete Ke; break; case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum: Ke=CreateKMatrixSlope(); - if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL); - delete Ke; break; case PrognosticAnalysisEnum: Ke=CreateKMatrixPrognostic(); - if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL); - delete Ke; break; case BalancedthicknessAnalysisEnum: Ke=CreateKMatrixBalancedthickness(); - if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL); - delete Ke; break; case BalancedvelocitiesAnalysisEnum: Ke=CreateKMatrixBalancedvelocities(); - if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL); - delete Ke; break; case ThermalAnalysisEnum: - CreateKMatrixThermal( Kgg); + Ke=CreateKMatrixThermal(); break; case MeltingAnalysisEnum: Ke=CreateKMatrixMelting(); - if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL); - delete Ke; 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; } } @@ -2840,25 +2830,31 @@ /*}}}*/ /*FUNCTION Penta::CreateKMatrixThermal {{{1*/ -void Penta::CreateKMatrixThermal(Mat Kgg){ - - /* local declarations */ +ElementMatrix* Penta::CreateKMatrixThermal(void){ + + /*compute all stiffness matrices for this element*/ + ElementMatrix* Ke1=CreateKMatrixThermalVolume(); + ElementMatrix* Ke2=CreateKMatrixThermalShelf(); + ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); + printf("-------------- file: Penta.cpp line: %i\n",__LINE__); + Ke->Echo(); + + /*clean-up and return*/ + delete Ke1; + delete Ke2; + return Ke; +} +/*}}}*/ +/*FUNCTION Penta::CreateKMatrixThermalVolume {{{1*/ +ElementMatrix* Penta::CreateKMatrixThermalVolume(void){ + int i,j; + int ig; int found=0; - - /* node data: */ const int numdof=NDOF1*NUMVERTICES; double xyz_list[NUMVERTICES][3]; int* doflist=NULL; - - /* gaussian points: */ - int ig; GaussPenta *gauss=NULL; - double K[2][2]={0.0}; - double u,v,w; - - /*matrices: */ - double K_terms[numdof][numdof]={0.0}; double Ke_gaussian_conduct[numdof][numdof]; double Ke_gaussian_advec[numdof][numdof]; @@ -2881,27 +2877,18 @@ double tBD_artdiff[3][numdof]; double tLD[numdof]; - double Jdet; - - /*Material properties: */ double gravity,rho_ice,rho_water; double heatcapacity,thermalconductivity; double mixed_layer_capacity,thermal_exchange_velocity; - - /*parameters: */ double dt,epsvel; bool artdiff; - - /*Collapsed formulation: */ Tria* tria=NULL; - /*If on water, skip: */ - if(IsOnWater())return; - - /* 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); - - // /*recovre material parameters: */ rho_water=matpar->GetRhoWater(); rho_ice=matpar->GetRhoIce(); @@ -2909,10 +2896,7 @@ heatcapacity=matpar->GetHeatCapacity(); thermalconductivity=matpar->GetThermalConductivity(); - - /*retrieve some parameters: */ this->parameters->FindParam(&dt,DtEnum); this->parameters->FindParam(&artdiff,ArtDiffEnum); this->parameters->FindParam(&epsvel,EpsVelEnum); - Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); @@ -2929,10 +2913,7 @@ /*Conduction: */ - /*Get B_conduct matrix: */ GetBConduct(&B_conduct[0][0],&xyz_list[0][0],gauss); - /*Build D: */ D_scalar=gauss->weight*Jdet*(thermalconductivity/(rho_ice*heatcapacity)); - if(dt) D_scalar=D_scalar*dt; @@ -2941,5 +2922,4 @@ D[2][0]=0; D[2][1]=0; D[2][2]=D_scalar; - /* Do the triple product B'*D*B: */ MatrixMultiply(&B_conduct[0][0],3,numdof,1,&D[0][0],3,3,0,&tBD_conduct[0][0],0); MatrixMultiply(&tBD_conduct[0][0],numdof,3,0,&B_conduct[0][0],3,numdof,0,&Ke_gaussian_conduct[0][0],0); @@ -2947,9 +2927,7 @@ /*Advection: */ - /*Get B_advec and Bprime_advec matrices: */ GetBAdvec(&B_advec[0][0],&xyz_list[0][0],gauss); GetBprimeAdvec(&Bprime_advec[0][0],&xyz_list[0][0],gauss); - //Build the D matrix vx_input->GetParameterValue(&u, gauss); vy_input->GetParameterValue(&v, gauss); @@ -2957,5 +2935,4 @@ D_scalar=gauss->weight*Jdet; - if(dt) D_scalar=D_scalar*dt; @@ -2964,9 +2941,9 @@ D[2][0]=0; D[2][1]=0; D[2][2]=D_scalar*w; - /* Do the triple product B'*D*Bprime: */ MatrixMultiply(&B_advec[0][0],3,numdof,1,&D[0][0],3,3,0,&tBD_advec[0][0],0); MatrixMultiply(&tBD_advec[0][0],numdof,3,0,&Bprime_advec[0][0],3,numdof,0,&Ke_gaussian_advec[0][0],0); /*Transient: */ + if(dt){ GetNodalFunctionsP1(&L[0], gauss); @@ -2974,5 +2951,4 @@ D_scalar=D_scalar; - /* Do the triple product L'*D*L: */ MatrixMultiply(&L[0],numdof,1,0,&D_scalar,1,1,0,&tLD[0],0); MatrixMultiply(&tLD[0],numdof,1,0,&L[0],1,numdof,0,&Ke_gaussian_transient[0][0],0); @@ -2983,4 +2959,5 @@ /*Artifficial diffusivity*/ + if(artdiff){ /*Build K: */ @@ -2990,8 +2967,6 @@ K[1][0]=D_scalar*fabs(u)*fabs(v);K[1][1]=D_scalar*pow(v,2); - /*Get B_artdiff: */ GetBArtdiff(&B_artdiff[0][0],&xyz_list[0][0],gauss); - /* Do the triple product B'*K*B: */ MatrixMultiply(&B_artdiff[0][0],2,numdof,1,&K[0][0],2,2,0,&tBD_artdiff[0][0],0); MatrixMultiply(&tBD_artdiff[0][0],numdof,2,0,&B_artdiff[0][0],2,numdof,0,&Ke_gaussian_artdiff[0][0],0); @@ -3001,23 +2976,23 @@ } - /*Add Ke_gaussian to pKe: */ - for(i=0;iCreateKMatrixThermal(Kgg); - delete tria->matice; delete tria; - } - - /*Free ressources:*/ + for(i=0;ivalues[i*numdof+j]+=Ke_gaussian_conduct[i][j]+Ke_gaussian_advec[i][j]+Ke_gaussian_transient[i][j]+Ke_gaussian_artdiff[i][j]; + } + + /*Clean up and return*/ delete gauss; - xfree((void**)&doflist); - + return Ke; +} +/*}}}*/ +/*FUNCTION Penta::CreateKMatrixThermalShelf {{{1*/ +ElementMatrix* Penta::CreateKMatrixThermalShelf(void){ + + if (!IsOnBed() || !IsOnShelf() || IsOnWater()) return NULL; + + /*Call Tria function*/ + Tria* tria=(Tria*)SpawnTria(0,1,2); + ElementMatrix* Ke=tria->CreateKMatrixThermal(); + delete tria->matice; delete tria; + + return Ke; } /*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 5878) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 5879) @@ -147,5 +147,7 @@ ElementMatrix* CreateKMatrixPrognostic(void); ElementMatrix* CreateKMatrixSlope(void); - void CreateKMatrixThermal(Mat Kggg); + ElementMatrix* CreateKMatrixThermal(void); + ElementMatrix* CreateKMatrixThermalVolume(void); + ElementMatrix* CreateKMatrixThermalShelf(void); void CreatePVectorBalancedthickness( Vec pg); void CreatePVectorBalancedvelocities( Vec pg); Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5878) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5879) @@ -3408,13 +3408,10 @@ /*}}}*/ /*FUNCTION Tria::CreateKMatrixThermal {{{1*/ -void Tria::CreateKMatrixThermal(Mat Kgg){ - +ElementMatrix* Tria::CreateKMatrixThermal(void){ + + 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; - double mixed_layer_capacity; double thermal_exchange_velocity; @@ -3422,9 +3419,5 @@ double rho_ice; double heatcapacity; - - int ig; GaussTria *gauss=NULL; - - /*matrices: */ double Jdet; double K_terms[numdof][numdof]={0.0}; @@ -3433,16 +3426,13 @@ double tl1l2l3D[3]; double D_scalar; - - /*parameters: */ double dt; - /*retrieve some parameters: */ + /*Initialize Element matrix and return if necessary*/ + if(IsOnWater() || !IsOnShelf()) return NULL; + ElementMatrix* Ke=this->NewElementMatrix(NoneApproximationEnum); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); this->parameters->FindParam(&dt,DtEnum); - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - - //recover material parameters mixed_layer_capacity=matpar->GetMixedLayerCapacity(); thermal_exchange_velocity=matpar->GetThermalExchangeVelocity(); @@ -3457,34 +3447,19 @@ gauss->GaussPoint(ig); - //Get the Jacobian determinant - GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0], gauss); - - /*Get nodal functions values: */ GetNodalFunctions(&l1l2l3[0], gauss); - /*Calculate DL on gauss point */ + GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0], gauss); D_scalar=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_velocity/(heatcapacity*rho_ice); - if(dt){ - D_scalar=dt*D_scalar; - } - - /* Do the triple product tL*D*L: */ + if(dt) D_scalar=dt*D_scalar; + MatrixMultiply(&l1l2l3[0],numdof,1,0,&D_scalar,1,1,0,&tl1l2l3D[0],0); MatrixMultiply(&tl1l2l3D[0],numdof,1,0,&l1l2l3[0],1,numdof,0,&Ke_gaussian[0][0],0); - for(i=0;i<3;i++){ - for(j=0;j<3;j++){ - K_terms[i][j]+=Ke_gaussian[i][j]; - } - } + for(i=0;ivalues[i*numdof+j]+=Ke_gaussian[i][j]; } - /*Add Ke_gg to global matrix Kgg: */ - MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)K_terms,ADD_VALUES); - - /*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 5878) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 5879) @@ -138,5 +138,5 @@ ElementMatrix* CreateKMatrixPrognostic_DG(void); ElementMatrix* CreateKMatrixSlope(void); - void CreateKMatrixThermal(Mat Kgg); + ElementMatrix* CreateKMatrixThermal(void); void CreatePVectorBalancedthickness_DG(Vec pg); void CreatePVectorBalancedthickness_CG(Vec pg);