Changeset 16007

Timestamp:

08/29/13 08:53:59 (12 years ago)

Author:

Mathieu Morlighem

Message:

NEW: added balancevelocity solver (not finished yet)

Location:

issm/trunk-jpl/src/c

Files:

• Makefile.am (modified) (1 diff)
• analyses/AnalysisConfiguration.cpp (modified) (1 diff)
• analyses/CorePointerFromSolutionEnum.cpp (modified) (1 diff)
• analyses/analyses.h (modified) (1 diff)
• analyses/balancevelocity_core.cpp (added)
• classes/Elements/Tria.cpp (modified) (6 diffs)
• classes/Elements/Tria.h (modified) (2 diffs)
• modules/ModelProcessorx/Balancevelocity (added)
• modules/ModelProcessorx/Balancevelocity/CreateConstraintsBalancevelocity.cpp (added)
• modules/ModelProcessorx/Balancevelocity/CreateLoadsBalancevelocity.cpp (added)
• modules/ModelProcessorx/Balancevelocity/CreateNodesBalancevelocity.cpp (added)
• modules/ModelProcessorx/Balancevelocity/UpdateElementsBalancevelocity.cpp (added)
• modules/ModelProcessorx/CreateDataSets.cpp (modified) (1 diff)
• modules/ModelProcessorx/DistributeNumDofs.cpp (modified) (1 diff)
• modules/ModelProcessorx/ModelProcessorx.h (modified) (1 diff)
• shared/Enum/EnumDefinitions.h (modified) (1 diff)
• shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• shared/Enum/StringToEnumx.cpp (modified) (3 diffs)

issm/trunk-jpl/src/c/Makefile.am ¶

@@ -501 +501 @@
                                             ./modules/ModelProcessorx/Balancethickness/CreateConstraintsBalancethickness.cpp\
                                                  ./modules/ModelProcessorx/Balancethickness/CreateLoadsBalancethickness.cpp\
+                                                 ./modules/ModelProcessorx/Balancevelocity/UpdateElementsBalancevelocity.cpp\
+                                                 ./modules/ModelProcessorx/Balancevelocity/CreateNodesBalancevelocity.cpp\
+                                                 ./modules/ModelProcessorx/Balancevelocity/CreateConstraintsBalancevelocity.cpp\
+                                                 ./modules/ModelProcessorx/Balancevelocity/CreateLoadsBalancevelocity.cpp\
                                                  ./analyses/balancethickness_core.cpp \
+                                                 ./analyses/balancevelocity_core.cpp \
                                                  ./analyses/dummy_core.cpp
 #}}}

issm/trunk-jpl/src/c/analyses/AnalysisConfiguration.cpp ¶

@@ -89 +89 @@
                         break;
 
+                case BalancevelocitySolutionEnum:
+                        numanalyses=2;
+                        analyses=xNew<int>(numanalyses);
+                        analyses[0]=BalancevelocityAnalysisEnum;
+                        analyses[1]=SurfaceSlopeAnalysisEnum;
+                        break;
+
                 case SurfaceSlopeSolutionEnum:
                         numanalyses=1;

issm/trunk-jpl/src/c/analyses/CorePointerFromSolutionEnum.cpp ¶

@@ -65 +65 @@
                         #endif
                         break;
+                case BalancevelocitySolutionEnum:
+                        #ifdef _HAVE_BALANCED_
+                        solutioncore=&balancevelocity_core;
+                        #else
+                        _error_("ISSM was not compiled with balanced capabilities. Exiting");
+                        #endif
+                        break;
                 case HydrologySolutionEnum:
                         #ifdef _HAVE_HYDROLOGY_

issm/trunk-jpl/src/c/analyses/analyses.h ¶

@@ -30 +30 @@
 void masstransport_core(FemModel* femmodel);
 void balancethickness_core(FemModel* femmodel);
+void balancevelocity_core(FemModel* femmodel);
 void slopecompute_core(FemModel* femmodel);
 void steadystate_core(FemModel* femmodel);

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp ¶

@@ -244 +244 @@
                         return CreateKMatrixBalancethickness();
                         break;
+                case BalancevelocityAnalysisEnum:
+                        return CreateKMatrixBalancevelocity();
+                        break;
                 #endif
                 #ifdef _HAVE_CONTROL_
@@ -385 +388 @@
                 case BalancethicknessAnalysisEnum:
                         return CreatePVectorBalancethickness();
+                        break;
+                case BalancevelocityAnalysisEnum:
+                        return CreatePVectorBalancevelocity();
                         break;
 #endif
@@ -1599 +1605 @@
                 case BalancethicknessAnalysisEnum:
                         InputUpdateFromSolutionOneDof(solution,ThicknessEnum);
+                        break;
+                case BalancevelocityAnalysisEnum:
+                        InputUpdateFromSolutionOneDof(solution,VelEnum);
                         break;
                 #endif
@@ -7067 +7076 @@
                 vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input);
         }
-        h=sqrt(2*this->GetArea());
+        h=sqrt(2.*this->GetArea());
 
         /*Start looping on the number of gaussian points:*/
@@ -7199 +7208 @@
 }
 /*}}}*/
+/*FUNCTION Tria::CreateKMatrixBalancevelocity{{{*/
+ElementMatrix* Tria::CreateKMatrixBalancevelocity(void){
+
+        /*Intermediaries */
+        IssmDouble xyz_list[NUMVERTICES][3];
+        IssmDouble dhdt_g,mb_g,ms_g,Jdettria;
+        IssmDouble h,gamma,thickness;
+        IssmDouble hnx,hny,dhnx[2],dhny[2];
+        int        i,j;
+        IssmDouble D_scalar;
+
+        /*Fetch number of nodes and dof 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<IssmDouble>(numnodes);
+        IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
+        IssmDouble*    HNx    = xNew<IssmDouble>(numnodes);
+        IssmDouble*    HNy    = xNew<IssmDouble>(numnodes);
+        IssmDouble*    H      = xNew<IssmDouble>(numnodes);
+        IssmDouble*    Nx     = xNew<IssmDouble>(numnodes);
+        IssmDouble*    Ny     = xNew<IssmDouble>(numnodes);
+
+        /*Retrieve all Inputs and parameters: */
+        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+        Input* H_input=inputs->GetInput(ThicknessEnum); _assert_(H_input);
+
+        /*Get vector N for all nodes*/
+        GetInputListOnNodes(Nx,SurfaceSlopeXEnum);
+        GetInputListOnNodes(Ny,SurfaceSlopeYEnum);
+        for(int i=0;i<numnodes;i++){
+                IssmDouble norm=sqrt(Nx[i]*Nx[i]+Ny[i]*Ny[i]);
+                Nx[i] = Nx[i]/norm;
+                Ny[i] = Ny[i]/norm;
+        }
+
+        /*Build HNx and HNy*/
+        GetInputListOnNodes(H, ThicknessEnum);
+        for(int i=0;i<numnodes;i++){
+                HNx[i]=H[i]*Nx[i];
+                HNy[i]=H[i]*Ny[i];
+        }
+
+        h=sqrt(2.*this->GetArea());
+
+        /*Start looping on the number of gaussian points:*/
+        GaussTria* gauss=new GaussTria(2);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+
+                H_input->GetInputValue(&thickness,gauss);
+                TriaRef::GetInputDerivativeValue(&dhnx[0],HNx,&xyz_list[0][0],gauss);
+                TriaRef::GetInputDerivativeValue(&dhny[0],HNy,&xyz_list[0][0],gauss);
+                TriaRef::GetInputValue(&hnx,HNx,gauss);
+                TriaRef::GetInputValue(&hny,HNy,gauss);
+
+                gamma=h/(2.*thickness);
+
+                GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
+                GetNodalFunctions(basis,gauss);
+                GetNodalFunctionsDerivatives(dbasis,&xyz_list[0][0],gauss);
+
+                for(int i=0;i<numnodes;i++){
+                        for(int j=0;j<numnodes;j++){
+                                Ke->values[i*numnodes+j] = gauss->weight*Jdettria*(
+                                                        basis[i]*basis[j]*(dhnx[0]+dhny[1])
+                                                        + basis[i]*(dbasis[0*numnodes+j]*hnx + dbasis[1*numnodes+j]*hny)
+                                                        + gamma*(basis[j]*(dhnx[0]+dhny[1]) + dbasis[0*numnodes+j]*hnx + dbasis[1*numnodes+j]*hny)*
+                                                        (basis[i]*(dhnx[0]+dhny[1]) + dbasis[0*numnodes+i]*hnx + dbasis[1*numnodes+i]*hny)
+                                                        );
+                        }
+                }
+        }
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(basis);
+        xDelete<IssmDouble>(dbasis);
+        xDelete<IssmDouble>(H);
+        xDelete<IssmDouble>(Nx);
+        xDelete<IssmDouble>(Ny);
+        xDelete<IssmDouble>(HNx);
+        xDelete<IssmDouble>(HNy);
+        delete gauss;
+        return Ke;
+}
+/*}}}*/
 /*FUNCTION Tria::CreatePVectorBalancethickness{{{*/
 ElementVector* Tria::CreatePVectorBalancethickness(void){
@@ -7293 +7390 @@
         /*Clean up and return*/
         xDelete<IssmDouble>(basis);
+        delete gauss;
+        return pe;
+}
+/*}}}*/
+/*FUNCTION Tria::CreatePVectorBalancevelocity{{{*/
+ElementVector* Tria::CreatePVectorBalancevelocity(void){
+
+        /*Intermediaries */
+        IssmDouble xyz_list[NUMVERTICES][3];
+        IssmDouble dhdt_g,mb_g,ms_g,Jdettria;
+        IssmDouble h,gamma,thickness;
+        IssmDouble hnx,hny,dhnx[2],dhny[2];
+
+        /*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<IssmDouble>(numnodes);
+        IssmDouble*    dbasis = xNew<IssmDouble>(numnodes*2);
+        IssmDouble*    HNx    = xNew<IssmDouble>(numnodes);
+        IssmDouble*    HNy    = xNew<IssmDouble>(numnodes);
+        IssmDouble*    H      = xNew<IssmDouble>(numnodes);
+        IssmDouble*    Nx     = xNew<IssmDouble>(numnodes);
+        IssmDouble*    Ny     = xNew<IssmDouble>(numnodes);
+
+        /*Retrieve all inputs and parameters*/
+        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+        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);
+        Input* H_input=inputs->GetInput(ThicknessEnum); _assert_(H_input);
+
+        /*Get vector N for all nodes*/
+        GetInputListOnNodes(Nx,SurfaceSlopeXEnum);
+        GetInputListOnNodes(Ny,SurfaceSlopeYEnum);
+        for(int i=0;i<numnodes;i++){
+                IssmDouble norm=sqrt(Nx[i]*Nx[i]+Ny[i]*Ny[i]);
+                Nx[i] = Nx[i]/norm;
+                Ny[i] = Ny[i]/norm;
+        }
+
+        /*Build HNx and HNy*/
+        GetInputListOnNodes(H, ThicknessEnum);
+        for(int i=0;i<numnodes;i++){
+                HNx[i]=H[i]*Nx[i];
+                HNy[i]=H[i]*Ny[i];
+        }
+
+        h=sqrt(2.*this->GetArea());
+
+        /* Start  looping on the number of gaussian points: */
+        GaussTria* gauss=new GaussTria(2);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+
+                ms_input->GetInputValue(&ms_g,gauss);
+                mb_input->GetInputValue(&mb_g,gauss);
+                dhdt_input->GetInputValue(&dhdt_g,gauss);
+                H_input->GetInputValue(&thickness,gauss);
+
+                GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
+                GetNodalFunctions(basis,gauss);
+                GetNodalFunctionsDerivatives(dbasis,&xyz_list[0][0],gauss);
+
+                TriaRef::GetInputDerivativeValue(&dhnx[0],HNx,&xyz_list[0][0],gauss);
+                TriaRef::GetInputDerivativeValue(&dhny[0],HNy,&xyz_list[0][0],gauss);
+                TriaRef::GetInputValue(&hnx,HNx,gauss);
+                TriaRef::GetInputValue(&hny,HNy,gauss);
+
+                gamma=h/(2.*thickness);
+
+                for(int i=0;i<numnodes;i++){
+                        pe->values[i]+=Jdettria*gauss->weight*(ms_g-mb_g-dhdt_g)*(
+                                                basis[i] + gamma*(basis[i]*(dhnx[0]+dhny[0])+hnx*dbasis[0] + hny*dbasis[1]));
+                }
+        }
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(basis);
+        xDelete<IssmDouble>(dbasis);
+        xDelete<IssmDouble>(H);
+        xDelete<IssmDouble>(Nx);
+        xDelete<IssmDouble>(Ny);
+        xDelete<IssmDouble>(HNx);
+        xDelete<IssmDouble>(HNy);
         delete gauss;
         return pe;

issm/trunk-jpl/src/c/classes/Elements/Tria.h ¶

@@ -183 +183 @@
                 ElementMatrix* CreateKMatrixBalancethickness_DG(void);
                 ElementMatrix* CreateKMatrixBalancethickness_CG(void);
+                ElementMatrix* CreateKMatrixBalancevelocity(void);
                 ElementMatrix* CreateKMatrixMelting(void);
                 ElementMatrix* CreateKMatrixMasstransport(void);
@@ -194 +195 @@
                 ElementVector* CreatePVectorBalancethickness_DG(void);
                 ElementVector* CreatePVectorBalancethickness_CG(void);
+                ElementVector* CreatePVectorBalancevelocity(void);
                 ElementVector* CreatePVectorMasstransport(void);
                 ElementVector* CreatePVectorMasstransport_CG(void);

issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateDataSets.cpp ¶

@@ -114 +114 @@
                         UpdateElementsBalancethickness(elements,iomodel,analysis_counter,analysis_type);
                         break;
+                case BalancevelocityAnalysisEnum:
+                        CreateNodesBalancevelocity(pnodes, iomodel);
+                        CreateConstraintsBalancevelocity(pconstraints,iomodel);
+                        CreateLoadsBalancevelocity(ploads,iomodel);
+                        UpdateElementsBalancevelocity(elements,iomodel,analysis_counter,analysis_type);
+                        break;
                 #endif
 

issm/trunk-jpl/src/c/modules/ModelProcessorx/DistributeNumDofs.cpp ¶

@@ -115 +115 @@
                         numdofs=1;
                         break;
+                case BalancevelocityAnalysisEnum:
+                        numdofs=1;
+                        break;
                 default:
                         _error_("analysis type: " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not implemented yet");

issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h ¶

@@ -119 +119 @@
 void UpdateElementsBalancethickness(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
 
+/*balancedvelocity:*/
+void CreateNodesBalancevelocity(Nodes** pnodes,IoModel* iomodel);
+void CreateConstraintsBalancevelocity(Constraints** pconstraints,IoModel* iomodel);
+void CreateLoadsBalancevelocity(Loads** ploads, IoModel* iomodel);
+void UpdateElementsBalancevelocity(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
+
 /*transient: */
 void UpdateElementsTransient(Elements* elements,Parameters* parameters,IoModel* iomodel,int analysis_counter,int analysis_type);

issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h ¶

@@ -273 +273 @@
         BalancethicknessSoftAnalysisEnum,
         BalancethicknessSoftSolutionEnum,
+        BalancevelocityAnalysisEnum,
+        BalancevelocitySolutionEnum,
         BedSlopeAnalysisEnum,
         BedSlopeSolutionEnum,

issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp ¶

@@ -279 +279 @@
                 case BalancethicknessSoftAnalysisEnum : return "BalancethicknessSoftAnalysis";
                 case BalancethicknessSoftSolutionEnum : return "BalancethicknessSoftSolution";
+                case BalancevelocityAnalysisEnum : return "BalancevelocityAnalysis";
+                case BalancevelocitySolutionEnum : return "BalancevelocitySolution";
                 case BedSlopeAnalysisEnum : return "BedSlopeAnalysis";
                 case BedSlopeSolutionEnum : return "BedSlopeSolution";

issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp ¶

@@ -285 +285 @@
               else if (strcmp(name,"BalancethicknessSoftAnalysis")==0) return BalancethicknessSoftAnalysisEnum;
               else if (strcmp(name,"BalancethicknessSoftSolution")==0) return BalancethicknessSoftSolutionEnum;
+              else if (strcmp(name,"BalancevelocityAnalysis")==0) return BalancevelocityAnalysisEnum;
+              else if (strcmp(name,"BalancevelocitySolution")==0) return BalancevelocitySolutionEnum;
               else if (strcmp(name,"BedSlopeAnalysis")==0) return BedSlopeAnalysisEnum;
               else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
@@ -381 +383 @@
               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 if (strcmp(name,"RiftfrontType")==0) return RiftfrontTypeEnum;
          else stage=4;
    }
    if(stage==4){
-              if (strcmp(name,"Segment")==0) return SegmentEnum;
+              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;
               else if (strcmp(name,"SpcDynamic")==0) return SpcDynamicEnum;
@@ -504 +506 @@
               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 if (strcmp(name,"P1P1")==0) return P1P1Enum;
          else stage=5;
    }
    if(stage==5){
-              if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum;
+              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;
               else if (strcmp(name,"MINIcondensed")==0) return MINIcondensedEnum;