Index: /issm/trunk-jpl/src/c/analyses/MeltingAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/MeltingAnalysis.cpp (revision 16839) +++ /issm/trunk-jpl/src/c/analyses/MeltingAnalysis.cpp (revision 16840) @@ -73,5 +73,46 @@ /*Finite Element Analysis*/ ElementMatrix* MeltingAnalysis::CreateKMatrix(Element* element){/*{{{*/ - _error_("not implemented yet"); + + /*Get basal element*/ + if(!element->IsOnBed()) return NULL; + Element* basalelement = element->SpawnBasalElement(); + + /*Intermediaries */ + IssmDouble D,Jdet; + IssmDouble *xyz_list = NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = basalelement->GetNumberOfNodes(); + + /*Initialize Element vector*/ + ElementMatrix* Ke = basalelement->NewElementMatrix(NoneApproximationEnum); + IssmDouble* basis = xNew(numnodes); + + /*Retrieve all inputs and parameters*/ + basalelement->GetVerticesCoordinates(&xyz_list); + IssmDouble latentheat = element->GetMaterialParameter(MaterialsLatentheatEnum); + IssmDouble heatcapacity = element->GetMaterialParameter(MaterialsHeatcapacityEnum); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=basalelement->NewGauss(2); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); + basalelement->NodalFunctions(basis,gauss); + D=latentheat/heatcapacity*gauss->weight*Jdet; + + TripleMultiply(basis,1,numnodes,1, + &D,1,1,0, + basis,1,numnodes,0, + &Ke->values[0],1); + } + + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(basis); + delete gauss; + basalelement->DeleteMaterials(); delete basalelement; + return Ke; }/*}}}*/ ElementVector* MeltingAnalysis::CreatePVector(Element* element){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp (revision 16839) +++ /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp (revision 16840) @@ -113,5 +113,188 @@ /*Finite Element Analysis*/ ElementMatrix* ThermalAnalysis::CreateKMatrix(Element* element){/*{{{*/ - _error_("not implemented yet"); + + /*compute all stiffness matrices for this element*/ + ElementMatrix* Ke1=CreateKMatrixVolume(element); + ElementMatrix* Ke2=CreateKMatrixShelf(element); + ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); + + /*clean-up and return*/ + delete Ke1; + delete Ke2; + return Ke; +}/*}}}*/ +ElementMatrix* ThermalAnalysis::CreateKMatrixVolume(Element* element){/*{{{*/ + + /*Intermediaries */ + int stabilization; + IssmDouble Jdet,dt,u,v,w,um,vm,wm,vel; + IssmDouble h,hx,hy,hz,vx,vy,vz; + IssmDouble tau_parameter,diameter; + IssmDouble D_scalar; + IssmDouble* xyz_list = NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Initialize Element vector and other vectors*/ + ElementMatrix* Ke = element->NewElementMatrix(); + IssmDouble* basis = xNew(numnodes); + IssmDouble* dbasis = xNew(3*numnodes); + IssmDouble* B = xNew(3*numnodes); + IssmDouble* Bprime = xNew(3*numnodes); + IssmDouble D[3][3]={0.}; + IssmDouble K[3][3]; + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + element->FindParam(&dt,TimesteppingTimeStepEnum); + element->FindParam(&stabilization,ThermalStabilizationEnum); + IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoWaterEnum); + IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); + IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum); + IssmDouble heatcapacity = element->GetMaterialParameter(MaterialsHeatcapacityEnum); + IssmDouble thermalconductivity = element->GetMaterialParameter(MaterialsThermalconductivityEnum); + IssmDouble kappa = thermalconductivity/(rho_ice*heatcapacity); + Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input = element->GetInput(VzEnum); _assert_(vz_input); + Input* vxm_input = element->GetInput(VxMeshEnum); _assert_(vxm_input); + Input* vym_input = element->GetInput(VyMeshEnum); _assert_(vym_input); + Input* vzm_input = element->GetInput(VzMeshEnum); _assert_(vzm_input); + if(stabilization==2) diameter=element->MinEdgeLength(xyz_list); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(2); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + D_scalar=gauss->weight*Jdet*kappa; + if(dt!=0.) D_scalar=D_scalar*dt; + + /*Conduction: */ + GetBConduct(B,element,xyz_list,gauss); + D[0][0]=D_scalar; + D[1][1]=D_scalar; + D[2][2]=D_scalar; + TripleMultiply(B,3,numnodes,1, + &D[0][0],3,3,0, + B,3,numnodes,0, + &Ke->values[0],1); + + /*Advection: */ + GetBAdvec(B,element,xyz_list,gauss); + GetBAdvecprime(Bprime,element,xyz_list,gauss); + vx_input->GetInputValue(&u,gauss); vxm_input->GetInputValue(&um,gauss); vx=u-um; + vy_input->GetInputValue(&v,gauss); vym_input->GetInputValue(&vm,gauss); vy=v-vm; + vz_input->GetInputValue(&w,gauss); vzm_input->GetInputValue(&wm,gauss); vz=w-wm; + D[0][0]=D_scalar*vx; + D[1][1]=D_scalar*vy; + D[2][2]=D_scalar*vz; + TripleMultiply(B,3,numnodes,1, + &D[0][0],3,3,0, + Bprime,3,numnodes,0, + &Ke->values[0],1); + + /*Transient: */ + if(dt!=0.){ + element->NodalFunctions(basis,gauss); + + TripleMultiply(basis,numnodes,1,0, + &D_scalar,1,1,0, + basis,1,numnodes,0, + &Ke->values[0],1); + } + + /*Artifficial diffusivity*/ + if(stabilization==1){ + element->ElementSizes(&hx,&hy,&hz); + vel=sqrt(vx*vx + vy*vy + vz*vz)+1.e-14; + h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) + pow(hz*vz/vel,2)); + K[0][0]=h/(2.*vel)*fabs(vx*vx); K[0][1]=h/(2.*vel)*fabs(vx*vy); K[0][2]=h/(2.*vel)*fabs(vx*vz); + K[1][0]=h/(2.*vel)*fabs(vy*vx); K[1][1]=h/(2.*vel)*fabs(vy*vy); K[1][2]=h/(2.*vel)*fabs(vy*vz); + K[2][0]=h/(2.*vel)*fabs(vz*vx); K[2][1]=h/(2.*vel)*fabs(vz*vy); K[2][2]=h/(2.*vel)*fabs(vz*vz); + for(int i=0;i<3;i++) for(int j=0;j<3;j++) K[i][j] = D_scalar*K[i][j]; + + GetBAdvecprime(Bprime,element,xyz_list,gauss); + + TripleMultiply(Bprime,3,numnodes,1, + &K[0][0],3,3,0, + Bprime,3,numnodes,0, + &Ke->values[0],1); + } + else if(stabilization==2){ + element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + tau_parameter=element->StabilizationParameter(u-um,v-vm,w-wm,diameter,kappa); + for(int i=0;ivalues[i*numnodes+j]+=tau_parameter*D_scalar* + ((u-um)*dbasis[0*3+i]+(v-vm)*dbasis[1*3+i]+(w-wm)*dbasis[2*3+i])*((u-um)*dbasis[0*3+j]+(v-vm)*dbasis[1*3+j]+(w-wm)*dbasis[2*3+j]); + } + } + if(dt!=0.){ + for(int i=0;ivalues[i*numnodes+j]+=tau_parameter*D_scalar*basis[j]*((u-um)*dbasis[0*3+i]+(v-vm)*dbasis[1*3+i]+(w-wm)*dbasis[2*3+i]); + } + } + } + } + } + + /*Clean up and return*/ + delete gauss; + return Ke; +}/*}}}*/ +ElementMatrix* ThermalAnalysis::CreateKMatrixShelf(Element* element){/*{{{*/ + + /*Initialize Element matrix and return if necessary*/ + if(!element->IsOnBed() || !element->IsFloating()) return NULL; + + IssmDouble dt,Jdet,D; + IssmDouble *xyz_list_base = NULL; + + /*Get basal element*/ + if(!element->IsOnBed() || !element->IsFloating()) return NULL; + + /*Fetch number of nodes for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Initialize vectors*/ + ElementMatrix* Ke = element->NewElementMatrix(); + IssmDouble* basis = xNew(numnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinatesBase(&xyz_list_base); + element->FindParam(&dt,TimesteppingTimeStepEnum); + IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum); + IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoWaterEnum); + IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); + IssmDouble heatcapacity = element->GetMaterialParameter(MaterialsHeatcapacityEnum); + IssmDouble mixed_layer_capacity= element->GetMaterialParameter(MaterialsMixedLayerCapacityEnum); + IssmDouble thermal_exchange_vel= element->GetMaterialParameter(MaterialsThermalExchangeVelocityEnum); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGaussBase(2); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); + element->NodalFunctions(basis,gauss); + + D=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_vel/(heatcapacity*rho_ice); + if(reCast(dt)) D=dt*D; + TripleMultiply(basis,numnodes,1,0, + &D,1,1,0, + basis,1,numnodes,0, + &Ke->values[0],1); + + } + + /*Clean up and return*/ + delete gauss; + xDelete(basis); + xDelete(xyz_list_base); + return Ke; }/*}}}*/ ElementVector* ThermalAnalysis::CreatePVector(Element* element){/*{{{*/ @@ -262,4 +445,91 @@ void ThermalAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ element->GetSolutionFromInputsOneDof(solution,TemperatureEnum); +}/*}}}*/ +void ThermalAnalysis::GetBConduct(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ + /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. + * For node i, Bi' can be expressed in the actual coordinate system + * by: + * Bi_conduct=[ dh/dx ] + * [ dh/dy ] + * [ dh/dz ] + * where h is the interpolation function for node i. + * + * We assume B has been allocated already, of size: 3x(NDOF1*numnodes) + */ + + /*Fetch number of nodes for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Get nodal functions derivatives*/ + IssmDouble* dbasis=xNew(3*numnodes); + element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + + /*Build B: */ + for(int i=0;i(dbasis); +}/*}}}*/ +void ThermalAnalysis::GetBAdvec(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ + /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. + * For node i, Bi' can be expressed in the actual coordinate system + * by: + * Bi_advec =[ h ] + * [ h ] + * [ h ] + * where h is the interpolation function for node i. + * + * We assume B has been allocated already, of size: 3x(NDOF1*NUMNODESP1) + */ + + /*Fetch number of nodes for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Get nodal functions*/ + IssmDouble* basis=xNew(numnodes); + element->NodalFunctions(basis,gauss); + + /*Build B: */ + for(int i=0;i(basis); +}/*}}}*/ +void ThermalAnalysis::GetBAdvecprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/ + /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. + * For node i, Bi' can be expressed in the actual coordinate system + * by: + * Biprime_advec=[ dh/dx ] + * [ dh/dy ] + * [ dh/dz ] + * where h is the interpolation function for node i. + * + * We assume B has been allocated already, of size: 3x(NDOF1*numnodes) + */ + + /*Fetch number of nodes for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Get nodal functions derivatives*/ + IssmDouble* dbasis=xNew(3*numnodes); + element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + + /*Build B: */ + for(int i=0;i(dbasis); }/*}}}*/ void ThermalAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.h (revision 16839) +++ /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.h (revision 16840) @@ -22,8 +22,13 @@ /*Finite element Analysis*/ ElementMatrix* CreateKMatrix(Element* element); + ElementMatrix* CreateKMatrixVolume(Element* element); + ElementMatrix* CreateKMatrixShelf(Element* element); ElementVector* CreatePVector(Element* element); ElementVector* CreatePVectorVolume(Element* element); ElementVector* CreatePVectorSheet(Element* element); ElementVector* CreatePVectorShelf(Element* element); + void GetBConduct(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); + void GetBAdvec(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); + void GetBAdvecprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); void GetSolutionFromInputs(Vector* solution,Element* element); void InputUpdateFromSolution(IssmDouble* solution,Element* element); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16839) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16840) @@ -49,4 +49,5 @@ virtual void CreateJacobianMatrix(Matrix* Jff)=0; virtual void DeleteMaterials(void)=0; + virtual void ElementSizes(IssmDouble* phx,IssmDouble* phy,IssmDouble* phz)=0; virtual void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure)=0; virtual IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16839) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16840) @@ -1210,6 +1210,6 @@ } /*}}}*/ -/*FUNCTION Penta::GetElementSizes{{{*/ -void Penta::GetElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){ +/*FUNCTION Penta::ElementSizes{{{*/ +void Penta::ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){ IssmDouble xyz_list[NUMVERTICES][3]; @@ -3981,5 +3981,5 @@ if(stabilization==1){ /*Build K: */ - GetElementSizes(&hx,&hy,&hz); + ElementSizes(&hx,&hy,&hz); vel=sqrt(vx*vx + vy*vy + vz*vz)+1.e-14; h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) + pow(hz*vz/vel,2)); @@ -4210,5 +4210,5 @@ if(stabilization==1){ /*Build K: */ - GetElementSizes(&hx,&hy,&hz); + ElementSizes(&hx,&hy,&hz); vel=sqrt(vx*vx + vy*vy + vz*vz)+1.e-14; h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) + pow(hz*vz/vel,2)); Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16839) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16840) @@ -75,4 +75,5 @@ void Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters); void DeleteMaterials(void){_error_("not implemented yet");}; + void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz); void FindParam(int* pvalue,int paramenum); void FindParam(IssmDouble* pvalue,int paramenum); @@ -221,5 +222,4 @@ IssmDouble GetGroundedPortion(IssmDouble* xyz_list); int GetElementType(void); - void GetElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz); Input* GetInput(int inputenum); void GetInputListOnVertices(IssmDouble* pvalue,int enumtype); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16839) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16840) @@ -84,4 +84,5 @@ void CreateJacobianMatrix(Matrix* Jff){_error_("not implemented yet");}; void Delta18oParameterization(void){_error_("not implemented yet");}; + void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");}; void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented yet");}; IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16839) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16840) @@ -79,4 +79,5 @@ void DeleteMaterials(void); void Delta18oParameterization(void); + void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");}; void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented yet");}; void FindParam(int* pvalue,int paramenum);