Index: /issm/trunk-jpl/src/c/analyses/balancevelocity_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/balancevelocity_core.cpp (revision 16009) +++ /issm/trunk-jpl/src/c/analyses/balancevelocity_core.cpp (revision 16010) @@ -20,5 +20,8 @@ if(VerboseSolution()) _printf0_("call computational core:\n"); - surfaceslope_core(femmodel); + femmodel->SetCurrentConfiguration(SurfaceSlopeAnalysisEnum,SmoothedSurfaceSlopeXAnalysisEnum); + solutionsequence_linear(femmodel); + femmodel->SetCurrentConfiguration(SurfaceSlopeAnalysisEnum,SmoothedSurfaceSlopeYAnalysisEnum); + solutionsequence_linear(femmodel); if(VerboseSolution()) _printf0_("call computational core:\n"); @@ -28,4 +31,6 @@ if(save_results){ if(VerboseSolution()) _printf0_(" saving results\n"); + InputToResultx(femmodel,SurfaceSlopeXEnum); + InputToResultx(femmodel,SurfaceSlopeYEnum); InputToResultx(femmodel,VelEnum); } Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16009) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16010) @@ -247,4 +247,7 @@ return CreateKMatrixBalancevelocity(); break; + case SmoothedSurfaceSlopeXAnalysisEnum: case SmoothedSurfaceSlopeYAnalysisEnum: + return CreateKMatrixSmoothedSlope(); + break; #endif #ifdef _HAVE_CONTROL_ @@ -391,4 +394,10 @@ case BalancevelocityAnalysisEnum: return CreatePVectorBalancevelocity(); + break; + case SmoothedSurfaceSlopeXAnalysisEnum: + return CreatePVectorSmoothedSlopeX(); + break; + case SmoothedSurfaceSlopeYAnalysisEnum: + return CreatePVectorSmoothedSlopeY(); break; #endif @@ -1608,4 +1617,10 @@ case BalancevelocityAnalysisEnum: InputUpdateFromSolutionOneDof(solution,VelEnum); + break; + case SmoothedSurfaceSlopeXAnalysisEnum: + InputUpdateFromSolutionOneDof(solution,SurfaceSlopeXEnum); + break; + case SmoothedSurfaceSlopeYAnalysisEnum: + InputUpdateFromSolutionOneDof(solution,SurfaceSlopeYEnum); break; #endif @@ -7296,4 +7311,52 @@ } /*}}}*/ +/*FUNCTION Tria::CreateKMatrixSmoothedSlope {{{*/ +ElementMatrix* Tria::CreateKMatrixSmoothedSlope(void){ + + /* Intermediaries */ + IssmDouble D_scalar,Jdet,thickness; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble D[2][2]; + IssmDouble l=12.; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = this->NumberofNodes(); + + /*Initialize Element matrix and vectors*/ + ElementMatrix* Ke = CreateMassMatrix(); + IssmDouble* B = xNew(2*numnodes); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); + + /* Start looping on the number of gaussian points: */ + GaussTria* gauss=new GaussTria(2); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetBMasstransport(B,&xyz_list[0][0], gauss); + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + thickness_input->GetInputValue(&thickness,gauss); + + D_scalar=gauss->weight*Jdet*(l*thickness)*(l*thickness); + D[0][0]=D_scalar; + D[1][0]=0.; + D[0][1]=0.; + D[1][1]=D_scalar; + + TripleMultiply(B,2,numnodes,1, + &D[0][0],2,2,0, + B,2,numnodes,0, + &Ke->values[0],1); + } + + /*Clean up and return*/ + delete gauss; + xDelete(B); + return Ke; +} +/*}}}*/ /*FUNCTION Tria::CreatePVectorBalancethickness{{{*/ ElementVector* Tria::CreatePVectorBalancethickness(void){ @@ -7477,4 +7540,90 @@ xDelete(HNx); xDelete(HNy); + delete gauss; + return pe; +} +/*}}}*/ +/*FUNCTION Tria::CreatePVectorSmoothedSlopeX{{{*/ +ElementVector* Tria::CreatePVectorSmoothedSlopeX(void){ + + /*Intermediaries */ + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble Jdettria; + IssmDouble thickness,slope[2]; + IssmDouble taud_x; + + /*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); + Input* H_input = inputs->GetInput(ThicknessEnum); _assert_(H_input); + Input* surface_input = inputs->GetInput(SurfaceEnum); _assert_(surface_input); + + /* Start looping on the number of gaussian points: */ + GaussTria* gauss=new GaussTria(2); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + H_input->GetInputValue(&thickness,gauss); + surface_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss); + taud_x = matpar->GetRhoIce()*matpar->GetG()*thickness*slope[0]; + + GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss); + GetNodalFunctions(basis,gauss); + + for(int i=0;ivalues[i]+=Jdettria*gauss->weight*taud_x*basis[i]; + } + + /*Clean up and return*/ + xDelete(basis); + delete gauss; + return pe; +} +/*}}}*/ +/*FUNCTION Tria::CreatePVectorSmoothedSlopeY{{{*/ +ElementVector* Tria::CreatePVectorSmoothedSlopeY(void){ + + /*Intermediaries */ + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble Jdettria; + IssmDouble thickness,slope[2]; + IssmDouble taud_y; + + /*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); + Input* H_input = inputs->GetInput(ThicknessEnum); _assert_(H_input); + Input* surface_input = inputs->GetInput(SurfaceEnum); _assert_(surface_input); + + /* Start looping on the number of gaussian points: */ + GaussTria* gauss=new GaussTria(2); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + H_input->GetInputValue(&thickness,gauss); + surface_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss); + taud_y = matpar->GetRhoIce()*matpar->GetG()*thickness*slope[1]; + + GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss); + GetNodalFunctions(basis,gauss); + + for(int i=0;ivalues[i]+=Jdettria*gauss->weight*taud_y*basis[i]; + } + + /*Clean up and return*/ + xDelete(basis); delete gauss; return pe; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16009) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16010) @@ -184,4 +184,5 @@ ElementMatrix* CreateKMatrixBalancethickness_CG(void); ElementMatrix* CreateKMatrixBalancevelocity(void); + ElementMatrix* CreateKMatrixSmoothedSlope(void); ElementMatrix* CreateKMatrixMelting(void); ElementMatrix* CreateKMatrixMasstransport(void); @@ -196,4 +197,6 @@ ElementVector* CreatePVectorBalancethickness_CG(void); ElementVector* CreatePVectorBalancevelocity(void); + ElementVector* CreatePVectorSmoothedSlopeX(void); + ElementVector* CreatePVectorSmoothedSlopeY(void); ElementVector* CreatePVectorMasstransport(void); ElementVector* CreatePVectorMasstransport_CG(void); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16009) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16010) @@ -301,4 +301,6 @@ SurfaceSlopeXAnalysisEnum, SurfaceSlopeYAnalysisEnum, + SmoothedSurfaceSlopeXAnalysisEnum, + SmoothedSurfaceSlopeYAnalysisEnum, ThermalAnalysisEnum, ThermalSolutionEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 16009) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 16010) @@ -307,4 +307,6 @@ case SurfaceSlopeXAnalysisEnum : return "SurfaceSlopeXAnalysis"; case SurfaceSlopeYAnalysisEnum : return "SurfaceSlopeYAnalysis"; + case SmoothedSurfaceSlopeXAnalysisEnum : return "SmoothedSurfaceSlopeXAnalysis"; + case SmoothedSurfaceSlopeYAnalysisEnum : return "SmoothedSurfaceSlopeYAnalysis"; case ThermalAnalysisEnum : return "ThermalAnalysis"; case ThermalSolutionEnum : return "ThermalSolution"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 16009) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 16010) @@ -313,4 +313,6 @@ else if (strcmp(name,"SurfaceSlopeXAnalysis")==0) return SurfaceSlopeXAnalysisEnum; else if (strcmp(name,"SurfaceSlopeYAnalysis")==0) return SurfaceSlopeYAnalysisEnum; + else if (strcmp(name,"SmoothedSurfaceSlopeXAnalysis")==0) return SmoothedSurfaceSlopeXAnalysisEnum; + else if (strcmp(name,"SmoothedSurfaceSlopeYAnalysis")==0) return SmoothedSurfaceSlopeYAnalysisEnum; else if (strcmp(name,"ThermalAnalysis")==0) return ThermalAnalysisEnum; else if (strcmp(name,"ThermalSolution")==0) return ThermalSolutionEnum; @@ -381,10 +383,10 @@ 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 stage=4; } if(stage==4){ - if (strcmp(name,"Riftfront")==0) return RiftfrontEnum; + 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; @@ -504,10 +506,10 @@ 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 stage=5; } if(stage==5){ - if (strcmp(name,"P1xP2")==0) return P1xP2Enum; + 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;