Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 8408) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 8409) @@ -531,4 +531,43 @@ }/*}}}*/ +/*FUNCTION Penta::CreateDVectorDiagnosticHoriz {{{1*/ +ElementVector* Penta::CreateDVectorDiagnosticHoriz(void){ + + int approximation; + inputs->GetParameterValue(&approximation,ApproximationEnum); + + switch(approximation){ + case StokesApproximationEnum: + return CreateDVectorDiagnosticStokes(); + default: + return NULL; //no need for doftypes outside of stokes approximation + } +} +/*}}}*/ +/*FUNCTION Penta::CreateDVectorDiagnosticStokes{{{1*/ +ElementVector* Penta::CreateDVectorDiagnosticStokes(void){ + + /*output: */ + ElementVector* De=NULL; + /*intermediary: */ + int approximation; + int i; + + /*Initialize Element vector and return if necessary*/ + inputs->GetParameterValue(&approximation,ApproximationEnum); + if(approximation!=StokesApproximationEnum) return NULL; + + De=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum); + + for (i=0;ivalues[i*4+0]=VelocityEnum; + De->values[i*4+1]=VelocityEnum; + De->values[i*4+2]=VelocityEnum; + De->values[i*4+3]=PressureEnum; + } + + return De; +} +/*}}}*/ /*FUNCTION Penta::CreateKMatrix {{{1*/ void Penta::CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs,Vec dg, Vec df){ @@ -1015,18 +1054,4 @@ } /*}}}*/ -/*FUNCTION Penta::CreateDVectorDiagnosticHoriz {{{1*/ -ElementVector* Penta::CreateDVectorDiagnosticHoriz(void){ - - int approximation; - inputs->GetParameterValue(&approximation,ApproximationEnum); - - switch(approximation){ - case StokesApproximationEnum: - return CreateDVectorDiagnosticStokes(); - default: - return NULL; //no need for doftypes outside of stokes approximation - } -} -/*}}}*/ /*FUNCTION Penta::CreateKMatrixDiagnosticHutter{{{1*/ ElementMatrix* Penta::CreateKMatrixDiagnosticHutter(void){ @@ -1388,29 +1413,4 @@ delete Ke2; return Ke; -} -/*}}}*/ -/*FUNCTION Penta::CreateDVectorDiagnosticStokes{{{1*/ -ElementVector* Penta::CreateDVectorDiagnosticStokes(void){ - - /*output: */ - ElementVector* De=NULL; - /*intermediary: */ - int approximation; - int i; - - /*Initialize Element vector and return if necessary*/ - inputs->GetParameterValue(&approximation,ApproximationEnum); - if(approximation!=StokesApproximationEnum) return NULL; - - De=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum); - - for (i=0;ivalues[i*4+0]=VelocityEnum; - De->values[i*4+1]=VelocityEnum; - De->values[i*4+2]=VelocityEnum; - De->values[i*4+3]=PressureEnum; - } - - return De; } /*}}}*/ @@ -1866,11 +1866,56 @@ ElementMatrix* Penta::CreateKMatrixThermalShelf(void){ + + /*Constants*/ + const int numdof=NDOF1*NUMVERTICES; + + /*Intermediaries */ + int i,j,ig; + double mixed_layer_capacity,thermal_exchange_velocity; + double rho_ice,rho_water,heatcapacity; + double Jdet2d,dt; + double xyz_list[NUMVERTICES][3]; + double xyz_list_tria[NUMVERTICES2D][3]; + double l1l6[NUMVERTICES]; + double D_scalar; + double Ke_gaussian[numdof][numdof]={0.0}; + GaussPenta *gauss=NULL; + + /*Initialize Element matrix and return if necessary*/ if (!IsOnBed() || !IsOnShelf()) return NULL; - - /*Call Tria function*/ - Tria* tria=(Tria*)SpawnTria(0,1,2); - ElementMatrix* Ke=tria->CreateKMatrixThermal(); - delete tria->matice; delete tria; - + ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum); + + /*Retrieve all inputs and parameters*/ + this->parameters->FindParam(&dt,DtEnum); + mixed_layer_capacity=matpar->GetMixedLayerCapacity(); + thermal_exchange_velocity=matpar->GetThermalExchangeVelocity(); + rho_water=matpar->GetRhoWater(); + rho_ice=matpar->GetRhoIce(); + heatcapacity=matpar->GetHeatCapacity(); + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + for(i=0;ibegin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss); + GetNodalFunctionsP1(&l1l6[0], gauss); + + D_scalar=gauss->weight*Jdet2d*rho_water*mixed_layer_capacity*thermal_exchange_velocity/(heatcapacity*rho_ice); + if(dt) D_scalar=dt*D_scalar; + + TripleMultiply(&l1l6[0],numdof,1,0, + &D_scalar,1,1,0, + &l1l6[0],1,numdof,0, + &Ke_gaussian[0][0],0); + + for(i=0;ivalues[i*numdof+j]+=Ke_gaussian[i][j]; + } + + /*Clean up and return*/ + delete gauss; return Ke; } @@ -2687,5 +2732,5 @@ water_pressure=gravity*rho_water*bed; - for(i=0;ivalues[i*NDOF4+j]+=(water_pressure+damper)*gauss->weight*Jdet2d*l1l6[i]*bed_normal[j]; + for(i=0;ivalues[i*NDOF4+j]+=(water_pressure+damper)*gauss->weight*Jdet2d*l1l6[i]*bed_normal[j]; } Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 8408) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 8409) @@ -1522,59 +1522,4 @@ } - /*Clean up and return*/ - delete gauss; - return Ke; -} -/*}}}*/ -/*FUNCTION Tria::CreateKMatrixThermal {{{1*/ -ElementMatrix* Tria::CreateKMatrixThermal(void){ - - /*Constants*/ - const int numdof=NDOF1*NUMVERTICES; - - /*Intermediaries */ - int i,j,ig; - double mixed_layer_capacity,thermal_exchange_velocity; - double rho_ice,rho_water,heatcapacity; - double Jdet,dt; - double xyz_list[NUMVERTICES][3]; - double l1l2l3[NUMVERTICES]; - double D_scalar; - double Ke_gaussian[numdof][numdof]={0.0}; - GaussTria *gauss=NULL; - - /*Initialize Element matrix and return if necessary*/ - if(!IsOnShelf()) return NULL; - ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum); - - /*Retrieve all inputs and parameters*/ - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - this->parameters->FindParam(&dt,DtEnum); - mixed_layer_capacity=matpar->GetMixedLayerCapacity(); - thermal_exchange_velocity=matpar->GetThermalExchangeVelocity(); - rho_water=matpar->GetRhoWater(); - rho_ice=matpar->GetRhoIce(); - heatcapacity=matpar->GetHeatCapacity(); - - /* Start looping on the number of gauss (nodes on the bedrock) */ - gauss=new GaussTria(2); - for (ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0], gauss); - GetNodalFunctions(&l1l2l3[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; - - TripleMultiply(&l1l2l3[0],numdof,1,0, - &D_scalar,1,1,0, - &l1l2l3[0],1,numdof,0, - &Ke_gaussian[0][0],0); - - for(i=0;ivalues[i*numdof+j]+=Ke_gaussian[i][j]; - } - /*Clean up and return*/ delete gauss; Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 8408) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 8409) @@ -160,5 +160,4 @@ ElementMatrix* CreateKMatrixPrognostic_DG(void); ElementMatrix* CreateKMatrixSlope(void); - ElementMatrix* CreateKMatrixThermal(void); ElementVector* CreatePVectorBalancethickness(void); ElementVector* CreatePVectorBalancethickness_DG(void);