Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 15876) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 15877) @@ -373,8 +373,17 @@ #Masstransport sources {{{ masstransport_sources = ./modules/ModelProcessorx/Masstransport/UpdateElementsMasstransport.cpp\ - ./modules/ModelProcessorx/Masstransport/CreateNodesMasstransport.cpp\ - ./modules/ModelProcessorx/Masstransport/CreateConstraintsMasstransport.cpp\ - ./modules/ModelProcessorx/Masstransport/CreateLoadsMasstransport.cpp\ - ./analyses/masstransport_core.cpp + ./modules/ModelProcessorx/Masstransport/CreateNodesMasstransport.cpp\ + ./modules/ModelProcessorx/Masstransport/CreateConstraintsMasstransport.cpp\ + ./modules/ModelProcessorx/Masstransport/CreateLoadsMasstransport.cpp\ + ./modules/ModelProcessorx/FreeSurfaceTop/UpdateElementsFreeSurfaceTop.cpp\ + ./modules/ModelProcessorx/FreeSurfaceTop/CreateNodesFreeSurfaceTop.cpp\ + ./modules/ModelProcessorx/FreeSurfaceTop/CreateConstraintsFreeSurfaceTop.cpp\ + ./modules/ModelProcessorx/FreeSurfaceTop/CreateLoadsFreeSurfaceTop.cpp\ + ./modules/ModelProcessorx/FreeSurfaceBase/UpdateElementsFreeSurfaceBase.cpp\ + ./modules/ModelProcessorx/FreeSurfaceBase/CreateNodesFreeSurfaceBase.cpp\ + ./modules/ModelProcessorx/FreeSurfaceBase/CreateConstraintsFreeSurfaceBase.cpp\ + ./modules/ModelProcessorx/FreeSurfaceBase/CreateLoadsFreeSurfaceBase.cpp\ + ./analyses/masstransport_core.cpp + #}}} #Thermal sources {{{ Index: /issm/trunk-jpl/src/c/analyses/AnalysisConfiguration.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/AnalysisConfiguration.cpp (revision 15876) +++ /issm/trunk-jpl/src/c/analyses/AnalysisConfiguration.cpp (revision 15877) @@ -108,15 +108,17 @@ case TransientSolutionEnum: - numanalyses=10-1; + numanalyses=11; analyses=xNew(numanalyses); - analyses[0]=StressbalanceAnalysisEnum; - analyses[1]=StressbalanceVerticalAnalysisEnum; - analyses[2]=StressbalanceSIAAnalysisEnum; - analyses[3]=SurfaceSlopeAnalysisEnum; - analyses[4]=BedSlopeAnalysisEnum; - analyses[5]=ThermalAnalysisEnum; - analyses[6]=MeltingAnalysisEnum; - analyses[7]=EnthalpyAnalysisEnum; - analyses[8]=MasstransportAnalysisEnum; + analyses[ 0]=StressbalanceAnalysisEnum; + analyses[ 1]=StressbalanceVerticalAnalysisEnum; + analyses[ 2]=StressbalanceSIAAnalysisEnum; + analyses[ 3]=SurfaceSlopeAnalysisEnum; + analyses[ 4]=BedSlopeAnalysisEnum; + analyses[ 5]=ThermalAnalysisEnum; + analyses[ 6]=MeltingAnalysisEnum; + analyses[ 7]=EnthalpyAnalysisEnum; + analyses[ 8]=MasstransportAnalysisEnum; + analyses[ 9]=FreeSurfaceBaseAnalysisEnum; + analyses[10]=FreeSurfaceTopAnalysisEnum; break; Index: /issm/trunk-jpl/src/c/analyses/masstransport_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/masstransport_core.cpp (revision 15876) +++ /issm/trunk-jpl/src/c/analyses/masstransport_core.cpp (revision 15877) @@ -13,11 +13,11 @@ /*parameters: */ - bool save_results; - bool issmbgradients,ispdd,isdelta18o; - int solution_type; + bool save_results; + bool issmbgradients,ispdd,isdelta18o,isFS,isfreesurface; + int solution_type; int *requested_outputs = NULL; - int numoutputs=0; + int numoutputs = 0; - /*activate formulation: */ + /*activate configuration*/ femmodel->SetCurrentConfiguration(MasstransportAnalysisEnum); @@ -27,4 +27,6 @@ femmodel->parameters->FindParam(&ispdd,SurfaceforcingsIspddEnum); femmodel->parameters->FindParam(&isdelta18o,SurfaceforcingsIsdelta18oEnum); + femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum); + femmodel->parameters->FindParam(&isfreesurface,MasstransportIsfreesurfaceEnum); femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); femmodel->parameters->FindParam(&numoutputs,MasstransportNumRequestedOutputsEnum); @@ -43,10 +45,22 @@ PositiveDegreeDayx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); } - if(VerboseSolution()) _printf0_(" call computational core\n"); - solutionsequence_linear(femmodel); + + if(isFS && isfreesurface){ + if(VerboseSolution()) _printf0_(" call free surface computational core\n"); + femmodel->SetCurrentConfiguration(FreeSurfaceBaseAnalysisEnum); + solutionsequence_linear(femmodel); + femmodel->SetCurrentConfiguration(FreeSurfaceTopAnalysisEnum); + solutionsequence_linear(femmodel); + } + else{ + if(VerboseSolution()) _printf0_(" call computational core\n"); + solutionsequence_linear(femmodel); + } if(save_results){ if(VerboseSolution()) _printf0_(" saving results\n"); InputToResultx(femmodel,ThicknessEnum); + InputToResultx(femmodel,BedEnum); + InputToResultx(femmodel,SurfaceEnum); femmodel->RequestedOutputsx(requested_outputs,numoutputs); } Index: /issm/trunk-jpl/src/c/analyses/transient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/transient_core.cpp (revision 15876) +++ /issm/trunk-jpl/src/c/analyses/transient_core.cpp (revision 15877) @@ -154,19 +154,14 @@ if(save_results){ if(VerboseSolution()) _printf0_(" saving transient results\n"); - InputToResultx(femmodel,SurfaceEnum); - InputToResultx(femmodel,BedEnum); InputToResultx(femmodel,SurfaceforcingsMassBalanceEnum); InputToResultx(femmodel,MaskElementonfloatingiceEnum); femmodel->RequestedOutputsx(requested_outputs,numoutputs); - if(isdelta18o){ InputToResultx(femmodel,SurfaceforcingsMonthlytemperaturesEnum); InputToResultx(femmodel,SurfaceforcingsPrecipitationEnum); - InputToResultx(femmodel,BasalFrictionEnum); } if(isgroundingline && (groundingline_migration==SubelementMigrationEnum || groundingline_migration==SubelementMigration2Enum)){ InputToResultx(femmodel,GLlevelsetEnum); } - if(VerboseSolution()) _printf0_(" saving temporary results\n"); OutputResultsx(femmodel->elements, femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,femmodel->results); Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15876) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15877) @@ -469,4 +469,10 @@ case MasstransportAnalysisEnum: return CreateKMatrixMasstransport(); + break; + case FreeSurfaceTopAnalysisEnum: + return CreateKMatrixFreeSurfaceTop(); + break; + case FreeSurfaceBaseAnalysisEnum: + return CreateKMatrixFreeSurfaceBase(); break; #ifdef _HAVE_BALANCED_ @@ -518,4 +524,30 @@ this->inputs->DeleteInput(VxAverageEnum); this->inputs->DeleteInput(VyAverageEnum); + + /*clean up and return*/ + return Ke; +} +/*}}}*/ +/*FUNCTION Penta::CreateKMatrixFreeSurfaceTop {{{*/ +ElementMatrix* Penta::CreateKMatrixFreeSurfaceTop(void){ + + if(!IsOnSurface()) return NULL; + + Tria* tria=(Tria*)SpawnTria(1); //lower face is 0, upper face is 1. + ElementMatrix* Ke=tria->CreateKMatrixFreeSurfaceTop(); + delete tria->material; delete tria; + + /*clean up and return*/ + return Ke; +} +/*}}}*/ +/*FUNCTION Penta::CreateKMatrixFreeSurfaceBase {{{*/ +ElementMatrix* Penta::CreateKMatrixFreeSurfaceBase(void){ + + if(!IsOnBed()) return NULL; + + Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1. + ElementMatrix* Ke=tria->CreateKMatrixFreeSurfaceBase(); + delete tria->material; delete tria; /*clean up and return*/ @@ -666,4 +698,10 @@ return CreatePVectorMasstransport(); break; + case FreeSurfaceTopAnalysisEnum: + return CreatePVectorFreeSurfaceTop(); + break; + case FreeSurfaceBaseAnalysisEnum: + return CreatePVectorFreeSurfaceBase(); + break; #ifdef _HAVE_BALANCED_ case BalancethicknessAnalysisEnum: @@ -704,4 +742,32 @@ this->inputs->DeleteInput(VxAverageEnum); this->inputs->DeleteInput(VyAverageEnum); + + /*Clean up and return*/ + return pe; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorFreeSurfaceTop {{{*/ +ElementVector* Penta::CreatePVectorFreeSurfaceTop(void){ + + if(!IsOnSurface()) return NULL; + + /*Call Tria function*/ + Tria* tria=(Tria*)SpawnTria(1); //lower face is 0, upper face is 1. + ElementVector* pe=tria->CreatePVectorFreeSurfaceTop(); + delete tria->material; delete tria; + + /*Clean up and return*/ + return pe; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorFreeSurfaceBase {{{*/ +ElementVector* Penta::CreatePVectorFreeSurfaceBase(void){ + + if(!IsOnBed()) return NULL; + + /*Call Tria function*/ + Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1. + ElementVector* pe=tria->CreatePVectorFreeSurfaceBase(); + delete tria->material; delete tria; /*Clean up and return*/ @@ -2263,4 +2329,10 @@ InputUpdateFromSolutionMasstransport(solution); break; + case FreeSurfaceTopAnalysisEnum: + InputUpdateFromSolutionFreeSurfaceTop(solution); + break; + case FreeSurfaceBaseAnalysisEnum: + InputUpdateFromSolutionFreeSurfaceBase(solution); + break; #ifdef _HAVE_BALANCED_ case BalancethicknessAnalysisEnum: @@ -2355,4 +2427,106 @@ /*Stop if we have reached the surface*/ if (penta->IsOnSurface()) break; + + /* get upper Penta*/ + penta=penta->GetUpperElement(); _assert_(penta->Id()!=this->id); + } + + /*Free ressources:*/ + xDelete(doflist); +} +/*}}}*/ +/*FUNCTION Penta::InputUpdateFromSolutionFreeSurfaceTop{{{*/ +void Penta::InputUpdateFromSolutionFreeSurfaceTop(IssmDouble* solution){ + + const int numdof = NDOF1*NUMVERTICES; + const int numdof2d = NDOF1*NUMVERTICES2D; + + int i; + int* doflist = NULL; + IssmDouble newthickness[numdof]; + IssmDouble newbed[numdof]; + IssmDouble newsurface[numdof]; + + /*If not on bed, return*/ + if (!IsOnSurface()) return; + + /*Get dof list: */ + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + + /*Use the dof list to index into the solution vector and extrude it */ + for(i=0;i(newsurface[i+numdof2d])) _error_("NaN found in solution vector"); + newsurface[i]=newsurface[i+numdof2d]; + } + + /*Get previous bed and thickness*/ + GetInputListOnVertices(&newbed[0],BedEnum); + + for(i=0;i0.); + } + + /*Start looping over all elements above current element and update all inputs*/ + Penta* penta=this; + for(;;){ + /*Add input to the element: */ + penta->inputs->AddInput(new PentaInput(ThicknessEnum,newthickness,P1Enum)); + penta->inputs->AddInput(new PentaInput(SurfaceEnum,newsurface,P1Enum)); + + /*Stop if we have reached the surface*/ + if(penta->IsOnBed()) break; + + /* get upper Penta*/ + penta=penta->GetLowerElement(); _assert_(penta->Id()!=this->id); + } + + /*Free ressources:*/ + xDelete(doflist); +} +/*}}}*/ +/*FUNCTION Penta::InputUpdateFromSolutionFreeSurfaceBase{{{*/ +void Penta::InputUpdateFromSolutionFreeSurfaceBase(IssmDouble* solution){ + + const int numdof = NDOF1*NUMVERTICES; + const int numdof2d = NDOF1*NUMVERTICES2D; + + int i; + int* doflist = NULL; + IssmDouble newthickness[numdof]; + IssmDouble newbed[numdof]; + IssmDouble newsurface[numdof]; + + /*If not on bed, return*/ + if (!IsOnBed()) return; + + /*Get dof list: */ + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + + /*Use the dof list to index into the solution vector and extrude it */ + for(i=0;i(newbed[i])) _error_("NaN found in solution vector"); + newbed[i+numdof2d]=newbed[i]; + } + + /*Get previous bed and thickness*/ + GetInputListOnVertices(&newsurface[0],SurfaceEnum); + + for(i=0;i0.); + } + + /*Start looping over all elements above current element and update all inputs*/ + Penta* penta=this; + for(;;){ + /*Add input to the element: */ + penta->inputs->AddInput(new PentaInput(ThicknessEnum,newthickness,P1Enum)); + penta->inputs->AddInput(new PentaInput(BedEnum,newbed,P1Enum)); + + /*Stop if we have reached the surface*/ + if(penta->IsOnSurface()) break; /* get upper Penta*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 15876) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 15877) @@ -181,6 +181,10 @@ ElementMatrix* CreateKMatrix(void); ElementMatrix* CreateKMatrixMasstransport(void); + ElementMatrix* CreateKMatrixFreeSurfaceTop(void); + ElementMatrix* CreateKMatrixFreeSurfaceBase(void); ElementVector* CreatePVector(void); ElementVector* CreatePVectorMasstransport(void); + ElementVector* CreatePVectorFreeSurfaceTop(void); + ElementVector* CreatePVectorFreeSurfaceBase(void); ElementVector* CreatePVectorSlope(void); void GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble xyz_zero[3][3],IssmDouble xyz_list[6][3],int numpoints); @@ -211,4 +215,6 @@ void InputExtrude(int enum_type,int object_type); void InputUpdateFromSolutionMasstransport(IssmDouble* solutiong); + void InputUpdateFromSolutionFreeSurfaceTop(IssmDouble* solutiong); + void InputUpdateFromSolutionFreeSurfaceBase(IssmDouble* solutiong); void InputUpdateFromSolutionOneDof(IssmDouble* solutiong,int enum_type); void InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solutiong,int enum_type); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15876) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15877) @@ -6487,4 +6487,210 @@ } /*}}}*/ +/*FUNCTION Tria::CreateKMatrixFreeSurfaceTop {{{*/ +ElementMatrix* Tria::CreateKMatrixFreeSurfaceTop(void){ + + /*Intermediaries */ + int stabilization; + int dim; + IssmDouble Jdettria,D_scalar,dt,h; + IssmDouble vel,vx,vy,dvxdx,dvydy; + IssmDouble dvx[2],dvy[2]; + IssmDouble xyz_list[NUMVERTICES][3]; + + /*Fetch number of nodes for this finite element*/ + int numnodes = this->NumberofNodes(); + + /*Initialize Element matrix and vectors*/ + ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum); + IssmDouble* basis = xNew(numnodes); + IssmDouble* B = xNew(2*numnodes); + IssmDouble* Bprime = xNew(2*numnodes); + IssmDouble D[2][2]; + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); + this->parameters->FindParam(&stabilization,MasstransportStabilizationEnum); + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + h=sqrt(2*this->GetArea()); + + /* Start looping on the number of gaussian points: */ + GaussTria *gauss=new GaussTria(2); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss); + GetNodalFunctions(basis,gauss); + + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + vx_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss); + vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss); + + D_scalar=gauss->weight*Jdettria; + + TripleMultiply(basis,1,numnodes,1, + &D_scalar,1,1,0, + basis,1,numnodes,0, + &Ke->values[0],1); + + GetBMasstransport(B,&xyz_list[0][0],gauss); + GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss); + + dvxdx=dvx[0]; + dvydy=dvy[1]; + D_scalar=dt*gauss->weight*Jdettria; + + D[0][0]=D_scalar*vx; + D[0][1]=0.; + D[1][0]=0.; + D[1][1]=D_scalar*vy; + TripleMultiply(B,2,numnodes,1, + &D[0][0],2,2,0, + Bprime,2,numnodes,0, + &Ke->values[0],1); + + if(stabilization==2){ + /*Streamline upwinding*/ + vel=sqrt(vx*vx+vy*vy)+1.e-8; + D[0][0]=h/(2*vel)*vx*vx; + D[1][0]=h/(2*vel)*vy*vx; + D[0][1]=h/(2*vel)*vx*vy; + D[1][1]=h/(2*vel)*vy*vy; + } + else if(stabilization==1){ + /*SSA*/ + vx_input->GetInputAverage(&vx); + vy_input->GetInputAverage(&vy); + D[0][0]=h/2.0*fabs(vx); + D[0][1]=0.; + D[1][0]=0.; + D[1][1]=h/2.0*fabs(vy); + } + if(stabilization==1 || stabilization==2){ + D[0][0]=D_scalar*D[0][0]; + D[1][0]=D_scalar*D[1][0]; + D[0][1]=D_scalar*D[0][1]; + D[1][1]=D_scalar*D[1][1]; + TripleMultiply(Bprime,2,numnodes,1, + &D[0][0],2,2,0, + Bprime,2,numnodes,0, + &Ke->values[0],1); + } + } + + /*Clean up and return*/ + xDelete(basis); + xDelete(B); + xDelete(Bprime); + delete gauss; + return Ke; +} +/*}}}*/ +/*FUNCTION Tria::CreateKMatrixFreeSurfaceBase {{{*/ +ElementMatrix* Tria::CreateKMatrixFreeSurfaceBase(void){ + + /*Intermediaries */ + int stabilization; + int dim; + IssmDouble Jdettria,D_scalar,dt,h; + IssmDouble vel,vx,vy,dvxdx,dvydy; + IssmDouble dvx[2],dvy[2]; + IssmDouble xyz_list[NUMVERTICES][3]; + + /*Fetch number of nodes for this finite element*/ + int numnodes = this->NumberofNodes(); + + /*Initialize Element matrix and vectors*/ + ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum); + IssmDouble* basis = xNew(numnodes); + IssmDouble* B = xNew(2*numnodes); + IssmDouble* Bprime = xNew(2*numnodes); + IssmDouble D[2][2]; + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); + this->parameters->FindParam(&stabilization,MasstransportStabilizationEnum); + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + h=sqrt(2*this->GetArea()); + + /* Start looping on the number of gaussian points: */ + GaussTria *gauss=new GaussTria(2); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss); + GetNodalFunctions(basis,gauss); + + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + vx_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss); + vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss); + + D_scalar=gauss->weight*Jdettria; + + TripleMultiply(basis,1,numnodes,1, + &D_scalar,1,1,0, + basis,1,numnodes,0, + &Ke->values[0],1); + + GetBMasstransport(B,&xyz_list[0][0],gauss); + GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss); + + dvxdx=dvx[0]; + dvydy=dvy[1]; + D_scalar=dt*gauss->weight*Jdettria; + + D[0][0]=D_scalar*vx; + D[0][1]=0.; + D[1][0]=0.; + D[1][1]=D_scalar*vy; + TripleMultiply(B,2,numnodes,1, + &D[0][0],2,2,0, + Bprime,2,numnodes,0, + &Ke->values[0],1); + + if(stabilization==2){ + /*Streamline upwinding*/ + vel=sqrt(vx*vx+vy*vy)+1.e-8; + D[0][0]=h/(2*vel)*vx*vx; + D[1][0]=h/(2*vel)*vy*vx; + D[0][1]=h/(2*vel)*vx*vy; + D[1][1]=h/(2*vel)*vy*vy; + } + else if(stabilization==1){ + /*SSA*/ + vx_input->GetInputAverage(&vx); + vy_input->GetInputAverage(&vy); + D[0][0]=h/2.0*fabs(vx); + D[0][1]=0.; + D[1][0]=0.; + D[1][1]=h/2.0*fabs(vy); + } + if(stabilization==1 || stabilization==2){ + D[0][0]=D_scalar*D[0][0]; + D[1][0]=D_scalar*D[1][0]; + D[0][1]=D_scalar*D[0][1]; + D[1][1]=D_scalar*D[1][1]; + TripleMultiply(Bprime,2,numnodes,1, + &D[0][0],2,2,0, + Bprime,2,numnodes,0, + &Ke->values[0],1); + } + } + + /*Clean up and return*/ + xDelete(basis); + xDelete(B); + xDelete(Bprime); + delete gauss; + return Ke; +} +/*}}}*/ /*FUNCTION Tria::CreatePVectorMasstransport{{{*/ ElementVector* Tria::CreatePVectorMasstransport(void){ @@ -6505,5 +6711,5 @@ /*Intermediaries */ IssmDouble Jdettria,dt; - IssmDouble surface_mass_balance_g,basal_melting_g,basal_melting_correction_g,thickness_g; + IssmDouble ms,mb,mb_correction,thickness; IssmDouble xyz_list[NUMVERTICES][3]; @@ -6518,10 +6724,10 @@ GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); - Input* surface_mass_balance_input = inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input); - Input* basal_melting_input = inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(basal_melting_input); - Input* basal_melting_correction_input = inputs->GetInput(BasalforcingsMeltingRateCorrectionEnum); - Input* thickness_input = inputs->GetInput(ThicknessEnum); _assert_(thickness_input); - - /*Initialize basal_melting_correction_g to 0, do not forget!:*/ + Input* ms_input = inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(ms_input); + Input* mb_input = inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(mb_input); + Input* mb_correction_input = inputs->GetInput(BasalforcingsMeltingRateCorrectionEnum); + Input* thickness_input = inputs->GetInput(ThicknessEnum); _assert_(thickness_input); + + /*Initialize mb_correction to 0, do not forget!:*/ /* Start looping on the number of gaussian points: */ GaussTria* gauss=new GaussTria(2); @@ -6533,13 +6739,13 @@ GetNodalFunctions(basis,gauss); - surface_mass_balance_input->GetInputValue(&surface_mass_balance_g,gauss); - basal_melting_input->GetInputValue(&basal_melting_g,gauss); - thickness_input->GetInputValue(&thickness_g,gauss); - if(basal_melting_correction_input) - basal_melting_correction_input->GetInputValue(&basal_melting_correction_g,gauss); + ms_input->GetInputValue(&ms,gauss); + mb_input->GetInputValue(&mb,gauss); + thickness_input->GetInputValue(&thickness,gauss); + if(mb_correction_input) + mb_correction_input->GetInputValue(&mb_correction,gauss); else - basal_melting_correction_g=0.; - - for(int i=0;ivalues[i]+=Jdettria*gauss->weight*(thickness_g+dt*(surface_mass_balance_g-basal_melting_g-basal_melting_correction_g))*basis[i]; + mb_correction=0.; + + for(int i=0;ivalues[i]+=Jdettria*gauss->weight*(thickness+dt*(ms-mb-mb_correction))*basis[i]; } @@ -6555,5 +6761,5 @@ /*Intermediaries */ IssmDouble Jdettria,dt; - IssmDouble surface_mass_balance_g,basal_melting_g,thickness_g; + IssmDouble ms,mb,thickness; IssmDouble xyz_list[NUMVERTICES][3]; @@ -6568,7 +6774,7 @@ this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* surface_mass_balance_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input); - Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(basal_melting_input); - Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* ms_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(ms_input); + Input* mb_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(mb_input); + Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); /* Start looping on the number of gaussian points: */ @@ -6581,9 +6787,102 @@ GetNodalFunctions(basis,gauss); - surface_mass_balance_input->GetInputValue(&surface_mass_balance_g,gauss); - basal_melting_input->GetInputValue(&basal_melting_g,gauss); - thickness_input->GetInputValue(&thickness_g,gauss); - - for(int i=0;ivalues[i]+=Jdettria*gauss->weight*(thickness_g+dt*(surface_mass_balance_g-basal_melting_g))*basis[i]; + ms_input->GetInputValue(&ms,gauss); + mb_input->GetInputValue(&mb,gauss); + thickness_input->GetInputValue(&thickness,gauss); + + for(int i=0;ivalues[i]+=Jdettria*gauss->weight*(thickness+dt*(ms-mb))*basis[i]; + } + + /*Clean up and return*/ + xDelete(basis); + delete gauss; + return pe; +} +/*}}}*/ +/*FUNCTION Tria::CreatePVectorFreeSurfaceTop {{{*/ +ElementVector* Tria::CreatePVectorFreeSurfaceTop(void){ + + /*Intermediaries */ + IssmDouble Jdettria,dt; + IssmDouble ms,surface,vz; + IssmDouble xyz_list[NUMVERTICES][3]; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = this->NumberofNodes(); + + /*Initialize Element vector and other vectors*/ + ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters); + IssmDouble* basis = xNew(numnodes); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); + Input* vz_input = inputs->GetInput(VzEnum); _assert_(vz_input); + Input* ms_input = inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(ms_input); + Input* surface_input= inputs->GetInput(ThicknessEnum); _assert_(surface_input); + + /*Initialize mb_correction to 0, do not forget!:*/ + /* Start looping on the number of gaussian points: */ + GaussTria* gauss=new GaussTria(2); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss); + GetNodalFunctions(basis,gauss); + + ms_input->GetInputValue(&ms,gauss); + surface_input->GetInputValue(&surface,gauss); + + for(int i=0;ivalues[i]+=Jdettria*gauss->weight*(surface + dt*ms + dt*vz)*basis[i]; + } + + /*Clean up and return*/ + xDelete(basis); + delete gauss; + return pe; +} +/*}}}*/ +/*FUNCTION Tria::CreatePVectorFreeSurfaceBase {{{*/ +ElementVector* Tria::CreatePVectorFreeSurfaceBase(void){ + + /*Intermediaries */ + IssmDouble Jdettria,dt; + IssmDouble mb,mb_correction,bed,vz; + IssmDouble xyz_list[NUMVERTICES][3]; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = this->NumberofNodes(); + + /*Initialize Element vector and other vectors*/ + ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters); + IssmDouble* basis = xNew(numnodes); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); + Input* vz_input = inputs->GetInput(VzEnum); _assert_(vz_input); + Input* mb_input = inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(mb_input); + Input* mb_correction_input = inputs->GetInput(BasalforcingsMeltingRateCorrectionEnum); + Input* bed_input = inputs->GetInput(BedEnum); _assert_(bed_input); + + /*Initialize mb_correction to 0, do not forget!:*/ + /* Start looping on the number of gaussian points: */ + GaussTria* gauss=new GaussTria(2); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss); + GetNodalFunctions(basis,gauss); + + mb_input->GetInputValue(&mb,gauss); + bed_input->GetInputValue(&bed,gauss); + if(mb_correction_input) + mb_correction_input->GetInputValue(&mb_correction,gauss); + else + mb_correction=0.; + + for(int i=0;ivalues[i]+=Jdettria*gauss->weight*(bed+dt*(mb-mb_correction) + dt*vz)*basis[i]; } @@ -6938,5 +7237,5 @@ /*Intermediaries */ IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble dhdt_g,basal_melting_g,surface_mass_balance_g,Jdettria; + IssmDouble dhdt_g,mb_g,ms_g,Jdettria; /*Fetch number of nodes and dof for this finite element*/ @@ -6949,7 +7248,7 @@ /*Retrieve all inputs and parameters*/ GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* surface_mass_balance_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input); - Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(basal_melting_input); - Input* dhdt_input=inputs->GetInput(BalancethicknessThickeningRateEnum); _assert_(dhdt_input); + Input* ms_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(ms_input); + Input* mb_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(mb_input); + Input* dhdt_input=inputs->GetInput(BalancethicknessThickeningRateEnum); _assert_(dhdt_input); /* Start looping on the number of gaussian points: */ @@ -6959,6 +7258,6 @@ gauss->GaussPoint(ig); - surface_mass_balance_input->GetInputValue(&surface_mass_balance_g,gauss); - basal_melting_input->GetInputValue(&basal_melting_g,gauss); + ms_input->GetInputValue(&ms_g,gauss); + mb_input->GetInputValue(&mb_g,gauss); dhdt_input->GetInputValue(&dhdt_g,gauss); @@ -6966,5 +7265,5 @@ GetNodalFunctions(basis,gauss); - for(int i=0;ivalues[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*basis[i]; + for(int i=0;ivalues[i]+=Jdettria*gauss->weight*(ms_g-mb_g-dhdt_g)*basis[i]; } @@ -6980,5 +7279,5 @@ /*Intermediaries */ IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble basal_melting_g,surface_mass_balance_g,dhdt_g,Jdettria; + IssmDouble mb_g,ms_g,dhdt_g,Jdettria; /*Fetch number of nodes and dof for this finite element*/ @@ -6991,7 +7290,7 @@ /*Retrieve all inputs and parameters*/ GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* surface_mass_balance_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input); - Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(basal_melting_input); - Input* dhdt_input=inputs->GetInput(BalancethicknessThickeningRateEnum); _assert_(dhdt_input); + Input* ms_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(ms_input); + Input* mb_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(mb_input); + Input* dhdt_input=inputs->GetInput(BalancethicknessThickeningRateEnum);_assert_(dhdt_input); /* Start looping on the number of gaussian points: */ @@ -7001,6 +7300,6 @@ gauss->GaussPoint(ig); - surface_mass_balance_input->GetInputValue(&surface_mass_balance_g,gauss); - basal_melting_input->GetInputValue(&basal_melting_g,gauss); + ms_input->GetInputValue(&ms_g,gauss); + mb_input->GetInputValue(&mb_g,gauss); dhdt_input->GetInputValue(&dhdt_g,gauss); @@ -7008,5 +7307,5 @@ GetNodalFunctions(basis,gauss); - for(int i=0;ivalues[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*basis[i]; + for(int i=0;ivalues[i]+=Jdettria*gauss->weight*(ms_g-mb_g-dhdt_g)*basis[i]; } Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 15876) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 15877) @@ -188,4 +188,6 @@ ElementMatrix* CreateKMatrixMasstransport_CG(void); ElementMatrix* CreateKMatrixMasstransport_DG(void); + ElementMatrix* CreateKMatrixFreeSurfaceTop(void); + ElementMatrix* CreateKMatrixFreeSurfaceBase(void); ElementMatrix* CreateMassMatrix(void); ElementVector* CreatePVector(void); @@ -196,4 +198,6 @@ ElementVector* CreatePVectorMasstransport_CG(void); ElementVector* CreatePVectorMasstransport_DG(void); + ElementVector* CreatePVectorFreeSurfaceTop(void); + ElementVector* CreatePVectorFreeSurfaceBase(void); ElementVector* CreatePVectorSlope(void); IssmDouble GetArea(void); Index: /issm/trunk-jpl/src/c/classes/Node.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Node.cpp (revision 15876) +++ /issm/trunk-jpl/src/c/classes/Node.cpp (revision 15877) @@ -93,4 +93,5 @@ /*2d solutions in 3d, we need to constrain all the nodes that are not on base*/ if( + analysis_type==FreeSurfaceBaseAnalysisEnum || analysis_type==MasstransportAnalysisEnum || analysis_type==MeltingAnalysisEnum || @@ -105,4 +106,13 @@ _assert_(iomodel->Data(MeshVertexonbedEnum)); if(!(reCast(iomodel->Data(MeshVertexonbedEnum)[io_index]))){ + this->Deactivate(); + } + } + } + if(analysis_type==FreeSurfaceTopAnalysisEnum){ + if(iomodel->dim==3){ + /*On a 3d mesh, we may have collapsed elements, hence dead nodes. Freeze them out: */ + _assert_(iomodel->Data(MeshVertexonsurfaceEnum)); + if(!(reCast(iomodel->Data(MeshVertexonsurfaceEnum)[io_index]))){ this->Deactivate(); } Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateDataSets.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateDataSets.cpp (revision 15876) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateDataSets.cpp (revision 15877) @@ -148,4 +148,16 @@ UpdateElementsMasstransport(elements,iomodel,analysis_counter,analysis_type); break; + case FreeSurfaceTopAnalysisEnum: + CreateNodesFreeSurfaceTop(pnodes, iomodel); + CreateConstraintsFreeSurfaceTop(pconstraints,iomodel); + CreateLoadsFreeSurfaceTop(ploads,iomodel); + UpdateElementsFreeSurfaceTop(elements,iomodel,analysis_counter,analysis_type); + break; + case FreeSurfaceBaseAnalysisEnum: + CreateNodesFreeSurfaceBase(pnodes, iomodel); + CreateConstraintsFreeSurfaceBase(pconstraints,iomodel); + CreateLoadsFreeSurfaceBase(ploads,iomodel); + UpdateElementsFreeSurfaceBase(elements,iomodel,analysis_counter,analysis_type); + break; #endif Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 15876) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 15877) @@ -57,4 +57,5 @@ parameters->AddObject(iomodel->CopyConstantObject(StressbalanceIsnewtonEnum)); parameters->AddObject(iomodel->CopyConstantObject(StressbalanceMaxiterEnum)); + parameters->AddObject(iomodel->CopyConstantObject(StressbalancePenaltyFactorEnum)); parameters->AddObject(iomodel->CopyConstantObject(SteadystateReltolEnum)); parameters->AddObject(iomodel->CopyConstantObject(SteadystateMaxiterEnum)); @@ -65,7 +66,7 @@ parameters->AddObject(iomodel->CopyConstantObject(TimesteppingTimeStepEnum)); parameters->AddObject(iomodel->CopyConstantObject(TimesteppingCflCoefficientEnum)); + parameters->AddObject(iomodel->CopyConstantObject(MasstransportIsfreesurfaceEnum)); parameters->AddObject(iomodel->CopyConstantObject(MasstransportHydrostaticAdjustmentEnum)); parameters->AddObject(iomodel->CopyConstantObject(MasstransportStabilizationEnum)); - parameters->AddObject(iomodel->CopyConstantObject(StressbalancePenaltyFactorEnum)); parameters->AddObject(iomodel->CopyConstantObject(MasstransportMinThicknessEnum)); parameters->AddObject(iomodel->CopyConstantObject(MasstransportPenaltyFactorEnum)); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/DistributeNumDofs.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/DistributeNumDofs.cpp (revision 15876) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/DistributeNumDofs.cpp (revision 15877) @@ -103,4 +103,10 @@ numdofs=1; break; + case FreeSurfaceTopAnalysisEnum: + numdofs=1; + break; + case FreeSurfaceBaseAnalysisEnum: + numdofs=1; + break; case GiaAnalysisEnum: numdofs=1; Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/FreeSurfaceBase/CreateConstraintsFreeSurfaceBase.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/FreeSurfaceBase/CreateConstraintsFreeSurfaceBase.cpp (revision 15877) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/FreeSurfaceBase/CreateConstraintsFreeSurfaceBase.cpp (revision 15877) @@ -0,0 +1,10 @@ +#include "../../../toolkits/toolkits.h" +#include "../../../classes/classes.h" +#include "../../../shared/shared.h" +#include "../ModelProcessorx.h" +#include "../../IoModelToConstraintsx/IoModelToConstraintsx.h" + +void CreateConstraintsFreeSurfaceBase(Constraints** pconstraints, IoModel* iomodel){ + + /*No constraints*/ +} Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/FreeSurfaceBase/CreateLoadsFreeSurfaceBase.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/FreeSurfaceBase/CreateLoadsFreeSurfaceBase.cpp (revision 15877) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/FreeSurfaceBase/CreateLoadsFreeSurfaceBase.cpp (revision 15877) @@ -0,0 +1,54 @@ +/*! \file CreateLoadsFreeSurfaceBase.c: + */ + +#include "../../../toolkits/toolkits.h" +#include "../../../classes/classes.h" +#include "../../../shared/shared.h" +#include "../ModelProcessorx.h" + +void CreateLoadsFreeSurfaceBase(Loads** ploads, IoModel* iomodel){ + + /*Intermediaries*/ + int element; + int penpair_ids[2]; + int count=0; + int numvertex_pairing; + + /*Recover pointer: */ + Loads* loads=*ploads; + + /*Create Penpair for vertex_pairing: */ + IssmDouble *vertex_pairing=NULL; + IssmDouble *nodeonbed=NULL; + iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,MasstransportVertexPairingEnum); + iomodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum); + for(int i=0;imy_vertices[reCast(vertex_pairing[2*i+0])-1]){ + + /*In debugging mode, check that the second node is in the same cpu*/ + _assert_(iomodel->my_vertices[reCast(vertex_pairing[2*i+1])-1]); + + /*Skip if one of the two is not on the bed*/ + if(!(reCast(nodeonbed[reCast(vertex_pairing[2*i+0])-1])) || !(reCast(nodeonbed[reCast(vertex_pairing[2*i+1])-1]))) continue; + + /*Get node ids*/ + penpair_ids[0]=iomodel->nodecounter+reCast(vertex_pairing[2*i+0]); + penpair_ids[1]=iomodel->nodecounter+reCast(vertex_pairing[2*i+1]); + + /*Create Load*/ + loads->AddObject(new Penpair( + iomodel->loadcounter+count+1, + &penpair_ids[0], + FreeSurfaceBaseAnalysisEnum)); + count++; + } + } + + /*free ressources: */ + iomodel->DeleteData(vertex_pairing,MasstransportVertexPairingEnum); + iomodel->DeleteData(nodeonbed,MeshVertexonbedEnum); + + /*Assign output pointer: */ + *ploads=loads; +} Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/FreeSurfaceBase/CreateNodesFreeSurfaceBase.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/FreeSurfaceBase/CreateNodesFreeSurfaceBase.cpp (revision 15877) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/FreeSurfaceBase/CreateNodesFreeSurfaceBase.cpp (revision 15877) @@ -0,0 +1,17 @@ +/* + * CreateNodesFreeSurfaceBase.c: + */ + +#include "../../../toolkits/toolkits.h" +#include "../../../classes/classes.h" +#include "../../../shared/shared.h" +#include "../../MeshPartitionx/MeshPartitionx.h" +#include "../ModelProcessorx.h" + +void CreateNodesFreeSurfaceBase(Nodes** pnodes, IoModel* iomodel){ + + /*Create Nodes either DG or CG depending on stabilization*/ + iomodel->FetchData(5,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,MaskVertexongroundediceEnum,MaskVertexonfloatingiceEnum,FlowequationVertexEquationEnum); + CreateNodes(pnodes,iomodel,FreeSurfaceBaseAnalysisEnum,P1Enum); + iomodel->DeleteData(5,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,MaskVertexongroundediceEnum,MaskVertexonfloatingiceEnum,FlowequationVertexEquationEnum); +} Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/FreeSurfaceBase/UpdateElementsFreeSurfaceBase.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/FreeSurfaceBase/UpdateElementsFreeSurfaceBase.cpp (revision 15877) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/FreeSurfaceBase/UpdateElementsFreeSurfaceBase.cpp (revision 15877) @@ -0,0 +1,41 @@ +/* + * UpdateElementsFreeSurfaceBase: + */ + +#include "../../../toolkits/toolkits.h" +#include "../../../classes/classes.h" +#include "../../../shared/shared.h" +#include "../../MeshPartitionx/MeshPartitionx.h" +#include "../ModelProcessorx.h" + +void UpdateElementsFreeSurfaceBase(Elements* elements, IoModel* iomodel,int analysis_counter,int analysis_type){ + + int finiteelement; + bool dakota_analysis; + + if(iomodel->dim!=3) _error_("not implemented yet for 2d models"); + /*Finite element type*/ + finiteelement = P1Enum; + + /*Update elements: */ + int counter=0; + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type,finiteelement); + counter++; + } + } + + iomodel->FetchDataToInput(elements,SurfaceEnum); + iomodel->FetchDataToInput(elements,MaskIcelevelsetEnum); + iomodel->FetchDataToInput(elements,MaskElementonfloatingiceEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum); + iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum); + iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateCorrectionEnum); + iomodel->FetchDataToInput(elements,VxEnum); + iomodel->FetchDataToInput(elements,VyEnum); + iomodel->FetchDataToInput(elements,VzEnum); + iomodel->FetchDataToInput(elements,MeshElementonbedEnum); + iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum); +} Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/FreeSurfaceTop/CreateConstraintsFreeSurfaceTop.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/FreeSurfaceTop/CreateConstraintsFreeSurfaceTop.cpp (revision 15877) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/FreeSurfaceTop/CreateConstraintsFreeSurfaceTop.cpp (revision 15877) @@ -0,0 +1,10 @@ +#include "../../../toolkits/toolkits.h" +#include "../../../classes/classes.h" +#include "../../../shared/shared.h" +#include "../ModelProcessorx.h" +#include "../../IoModelToConstraintsx/IoModelToConstraintsx.h" + +void CreateConstraintsFreeSurfaceTop(Constraints** pconstraints, IoModel* iomodel){ + + /*No constraints*/ +} Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/FreeSurfaceTop/CreateLoadsFreeSurfaceTop.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/FreeSurfaceTop/CreateLoadsFreeSurfaceTop.cpp (revision 15877) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/FreeSurfaceTop/CreateLoadsFreeSurfaceTop.cpp (revision 15877) @@ -0,0 +1,54 @@ +/*! \file CreateLoadsFreeSurfaceTop.c: + */ + +#include "../../../toolkits/toolkits.h" +#include "../../../classes/classes.h" +#include "../../../shared/shared.h" +#include "../ModelProcessorx.h" + +void CreateLoadsFreeSurfaceTop(Loads** ploads, IoModel* iomodel){ + + /*Intermediaries*/ + int element; + int penpair_ids[2]; + int count=0; + int numvertex_pairing; + + /*Recover pointer: */ + Loads* loads=*ploads; + + /*Create Penpair for vertex_pairing: */ + IssmDouble *vertex_pairing=NULL; + IssmDouble *nodeonsurface=NULL; + iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,MasstransportVertexPairingEnum); + iomodel->FetchData(&nodeonsurface,NULL,NULL,MeshVertexonsurfaceEnum); + for(int i=0;imy_vertices[reCast(vertex_pairing[2*i+0])-1]){ + + /*In debugging mode, check that the second node is in the same cpu*/ + _assert_(iomodel->my_vertices[reCast(vertex_pairing[2*i+1])-1]); + + /*Skip if one of the two is not on the bed*/ + if(!(reCast(nodeonsurface[reCast(vertex_pairing[2*i+0])-1])) || !(reCast(nodeonsurface[reCast(vertex_pairing[2*i+1])-1]))) continue; + + /*Get node ids*/ + penpair_ids[0]=iomodel->nodecounter+reCast(vertex_pairing[2*i+0]); + penpair_ids[1]=iomodel->nodecounter+reCast(vertex_pairing[2*i+1]); + + /*Create Load*/ + loads->AddObject(new Penpair( + iomodel->loadcounter+count+1, + &penpair_ids[0], + FreeSurfaceTopAnalysisEnum)); + count++; + } + } + + /*free ressources: */ + iomodel->DeleteData(vertex_pairing,MasstransportVertexPairingEnum); + iomodel->DeleteData(nodeonsurface,MeshVertexonsurfaceEnum); + + /*Assign output pointer: */ + *ploads=loads; +} Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/FreeSurfaceTop/CreateNodesFreeSurfaceTop.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/FreeSurfaceTop/CreateNodesFreeSurfaceTop.cpp (revision 15877) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/FreeSurfaceTop/CreateNodesFreeSurfaceTop.cpp (revision 15877) @@ -0,0 +1,17 @@ +/* + * CreateNodesFreeSurfaceTop.c: + */ + +#include "../../../toolkits/toolkits.h" +#include "../../../classes/classes.h" +#include "../../../shared/shared.h" +#include "../../MeshPartitionx/MeshPartitionx.h" +#include "../ModelProcessorx.h" + +void CreateNodesFreeSurfaceTop(Nodes** pnodes, IoModel* iomodel){ + + /*Create Nodes either DG or CG depending on stabilization*/ + iomodel->FetchData(5,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,MaskVertexongroundediceEnum,MaskVertexonfloatingiceEnum,FlowequationVertexEquationEnum); + CreateNodes(pnodes,iomodel,FreeSurfaceTopAnalysisEnum,P1Enum); + iomodel->DeleteData(5,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,MaskVertexongroundediceEnum,MaskVertexonfloatingiceEnum,FlowequationVertexEquationEnum); +} Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/FreeSurfaceTop/UpdateElementsFreeSurfaceTop.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/FreeSurfaceTop/UpdateElementsFreeSurfaceTop.cpp (revision 15877) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/FreeSurfaceTop/UpdateElementsFreeSurfaceTop.cpp (revision 15877) @@ -0,0 +1,41 @@ +/* + * UpdateElementsFreeSurfaceTop: + */ + +#include "../../../toolkits/toolkits.h" +#include "../../../classes/classes.h" +#include "../../../shared/shared.h" +#include "../../MeshPartitionx/MeshPartitionx.h" +#include "../ModelProcessorx.h" + +void UpdateElementsFreeSurfaceTop(Elements* elements, IoModel* iomodel,int analysis_counter,int analysis_type){ + + int finiteelement; + bool dakota_analysis; + + if(iomodel->dim!=3) _error_("not implemented yet for 2d models"); + /*Finite element type*/ + finiteelement = P1Enum; + + /*Update elements: */ + int counter=0; + for(int i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + Element* element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type,finiteelement); + counter++; + } + } + + iomodel->FetchDataToInput(elements,SurfaceEnum); + iomodel->FetchDataToInput(elements,MaskIcelevelsetEnum); + iomodel->FetchDataToInput(elements,MaskElementonfloatingiceEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum); + iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum); + iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateCorrectionEnum); + iomodel->FetchDataToInput(elements,VxEnum); + iomodel->FetchDataToInput(elements,VyEnum); + iomodel->FetchDataToInput(elements,VzEnum); + iomodel->FetchDataToInput(elements,MeshElementonbedEnum); + iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum); +} Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h (revision 15876) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h (revision 15877) @@ -104,4 +104,12 @@ void CreateLoadsMasstransport(Loads** ploads, IoModel* iomodel); void UpdateElementsMasstransport(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); +void CreateNodesFreeSurfaceTop(Nodes** pnodes,IoModel* iomodel); +void CreateConstraintsFreeSurfaceTop(Constraints** pconstraints,IoModel* iomodel); +void CreateLoadsFreeSurfaceTop(Loads** ploads, IoModel* iomodel); +void UpdateElementsFreeSurfaceTop(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); +void CreateNodesFreeSurfaceBase(Nodes** pnodes,IoModel* iomodel); +void CreateConstraintsFreeSurfaceBase(Constraints** pconstraints,IoModel* iomodel); +void CreateLoadsFreeSurfaceBase(Loads** ploads, IoModel* iomodel); +void UpdateElementsFreeSurfaceBase(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); /*balancedthickness:*/ Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 15876) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 15877) @@ -188,4 +188,5 @@ MiscellaneousNameEnum, //FIXME: only used by qmu, should not be marshalled (already in queueing script) MasstransportHydrostaticAdjustmentEnum, + MasstransportIsfreesurfaceEnum, MasstransportMinThicknessEnum, MasstransportPenaltyFactorEnum, @@ -295,4 +296,6 @@ MasstransportAnalysisEnum, MasstransportSolutionEnum, + FreeSurfaceBaseAnalysisEnum, + FreeSurfaceTopAnalysisEnum, SteadystateSolutionEnum, SurfaceSlopeAnalysisEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 15876) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 15877) @@ -196,4 +196,5 @@ case MiscellaneousNameEnum : return "MiscellaneousName"; case MasstransportHydrostaticAdjustmentEnum : return "MasstransportHydrostaticAdjustment"; + case MasstransportIsfreesurfaceEnum : return "MasstransportIsfreesurface"; case MasstransportMinThicknessEnum : return "MasstransportMinThickness"; case MasstransportPenaltyFactorEnum : return "MasstransportPenaltyFactor"; @@ -301,4 +302,6 @@ case MasstransportAnalysisEnum : return "MasstransportAnalysis"; case MasstransportSolutionEnum : return "MasstransportSolution"; + case FreeSurfaceBaseAnalysisEnum : return "FreeSurfaceBaseAnalysis"; + case FreeSurfaceTopAnalysisEnum : return "FreeSurfaceTopAnalysis"; case SteadystateSolutionEnum : return "SteadystateSolution"; case SurfaceSlopeAnalysisEnum : return "SurfaceSlopeAnalysis"; @@ -587,4 +590,5 @@ case PatersonEnum : return "Paterson"; case ArrheniusEnum : return "Arrhenius"; + case LliboutryDuvalEnum : return "LliboutryDuval"; case MaximumNumberOfDefinitionsEnum : return "MaximumNumberOfDefinitions"; default : return "unknown"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 15876) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 15877) @@ -199,4 +199,5 @@ else if (strcmp(name,"MiscellaneousName")==0) return MiscellaneousNameEnum; else if (strcmp(name,"MasstransportHydrostaticAdjustment")==0) return MasstransportHydrostaticAdjustmentEnum; + else if (strcmp(name,"MasstransportIsfreesurface")==0) return MasstransportIsfreesurfaceEnum; else if (strcmp(name,"MasstransportMinThickness")==0) return MasstransportMinThicknessEnum; else if (strcmp(name,"MasstransportPenaltyFactor")==0) return MasstransportPenaltyFactorEnum; @@ -259,9 +260,9 @@ else if (strcmp(name,"GiaMantleViscosity")==0) return GiaMantleViscosityEnum; else if (strcmp(name,"GiaLithosphereThickness")==0) return GiaLithosphereThicknessEnum; - else if (strcmp(name,"Thickness")==0) return ThicknessEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"TimesteppingStartTime")==0) return TimesteppingStartTimeEnum; + if (strcmp(name,"Thickness")==0) return ThicknessEnum; + else if (strcmp(name,"TimesteppingStartTime")==0) return TimesteppingStartTimeEnum; else if (strcmp(name,"TimesteppingFinalTime")==0) return TimesteppingFinalTimeEnum; else if (strcmp(name,"TimesteppingCflCoefficient")==0) return TimesteppingCflCoefficientEnum; @@ -307,4 +308,6 @@ else if (strcmp(name,"MasstransportAnalysis")==0) return MasstransportAnalysisEnum; else if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum; + else if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum; + else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum; else if (strcmp(name,"SteadystateSolution")==0) return SteadystateSolutionEnum; else if (strcmp(name,"SurfaceSlopeAnalysis")==0) return SurfaceSlopeAnalysisEnum; @@ -380,11 +383,11 @@ else if (strcmp(name,"Profiler")==0) return ProfilerEnum; else if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum; - else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum; - else if (strcmp(name,"VectorParam")==0) return VectorParamEnum; - else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"RiftfrontType")==0) return RiftfrontTypeEnum; + if (strcmp(name,"NodeSId")==0) return NodeSIdEnum; + else if (strcmp(name,"VectorParam")==0) return VectorParamEnum; + else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum; + else if (strcmp(name,"RiftfrontType")==0) return RiftfrontTypeEnum; else if (strcmp(name,"Segment")==0) return SegmentEnum; else if (strcmp(name,"SegmentRiftfront")==0) return SegmentRiftfrontEnum; @@ -503,11 +506,11 @@ else if (strcmp(name,"P1bubble")==0) return P1bubbleEnum; else if (strcmp(name,"P1bubblecondensed")==0) return P1bubblecondensedEnum; - else if (strcmp(name,"P2")==0) return P2Enum; - else if (strcmp(name,"P2xP1")==0) return P2xP1Enum; - else if (strcmp(name,"P1xP2")==0) return P1xP2Enum; else stage=5; } if(stage==5){ - if (strcmp(name,"P1P1")==0) return P1P1Enum; + if (strcmp(name,"P2")==0) return P2Enum; + else if (strcmp(name,"P2xP1")==0) return P2xP1Enum; + else if (strcmp(name,"P1xP2")==0) return P1xP2Enum; + else if (strcmp(name,"P1P1")==0) return P1P1Enum; else if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum; else if (strcmp(name,"MINI")==0) return MINIEnum; @@ -599,4 +602,5 @@ else if (strcmp(name,"Paterson")==0) return PatersonEnum; else if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum; + else if (strcmp(name,"LliboutryDuval")==0) return LliboutryDuvalEnum; else if (strcmp(name,"MaximumNumberOfDefinitions")==0) return MaximumNumberOfDefinitionsEnum; else stage=6;