source: issm/oecreview/Archive/17984-18295/ISSM-18283-18284.diff@ 18296

../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h

@@ -78 +78 @@
                 ElementVector* CreatePVectorFSViscousXTH(Element* element);
                 ElementVector* CreatePVectorFSShelf(Element* element);
                 ElementVector* CreatePVectorFSFront(Element* element);
+                ElementVector* CreatePVectorFSFriction(Element* element);
+                ElementVector* CreatePVectorFSStress(Element* element);
                 void GetBFS(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
                 void GetBFSprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
                 void GetBFSvel(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);

../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

@@ -244 +244 @@
         iomodel->Constant(&fe_FS,FlowequationFeFSEnum);
         if(fe_FS==LATaylorHoodEnum){
                 iomodel->FetchDataToInput(elements,PressureEnum);
+                InputUpdateFromConstantx(elements,0.,SigmaNNEnum);
         }
 
         /*Friction law variables*/
@@ -2895 +2896 @@
 
         /*Fetch number of nodes and dof for this finite element*/
         int vnumnodes = element->NumberofNodesVelocity();
-        int pnumnodes;
-        if(dim==2) pnumnodes=3;
-        else pnumnodes=6;
-        int lnumnodes;
-        if(dim==2) lnumnodes=2;
-        else lnumnodes=3;
-        //int pnumnodes = element->NumberofNodes(P1Enum);
+        int pnumnodes = element->GetNumberOfNodes(P1Enum);
+        int lnumnodes = element->GetNumberOfNodes(P2Enum);
         int numdof    = vnumnodes*dim;
         int pnumdof   = pnumnodes;
         int lnumdof   = lnumnodes;
@@ -2919 +2915 @@
         IssmDouble*    BU       = xNew<IssmDouble>(pnumdof);
         IssmDouble*    BprimeU  = xNew<IssmDouble>(numdof);
         IssmDouble*    D        = xNewZeroInit<IssmDouble>(epssize*epssize);
-        IssmDouble*    CtCUzawa = xNewZeroInit<IssmDouble>(numdof*pnumdof);
-        IssmDouble*    C        = xNew<IssmDouble>(lnumdof);
-        IssmDouble*    Cprime   = xNew<IssmDouble>(numdof);
 
         /*Retrieve all inputs and parameters*/
         element->GetVerticesCoordinates(&xyz_list);
-        Input* vx_input=element->GetInput(VxEnum);     _assert_(vx_input);
-        Input* vy_input=element->GetInput(VyEnum);     _assert_(vy_input);
-        Input* vz_input;
-        if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);}
+        Input* vx_input = element->GetInput(VxEnum);     _assert_(vx_input);
+        Input* vy_input = element->GetInput(VyEnum);     _assert_(vy_input);
+        Input* vz_input = NULL;
+        if(dim==3){vz_input = element->GetInput(VzEnum); _assert_(vz_input);}
 
         /* Start  looping on the number of gaussian points: */
         Gauss* gauss = element->NewGauss(5);
@@ -2956 +2949 @@
                                         &DU,1,1,0,
                                         BprimeU,1,numdof,0,
                                         BtBUzawa,1);
-
         }
 
         if(element->IsOnBase()){ 
@@ -2964 +2956 @@
                 element->GetVerticesCoordinatesBase(&xyz_list_base);
                 element->NormalBase(&normal[0],xyz_list_base);
 
-                int         lsize;
-                IssmDouble* Dlambda = xNewZeroInit<IssmDouble>(dim*dim);
-                IssmDouble* C       = xNewZeroInit<IssmDouble>(dim*lnumdof);
-                IssmDouble* Cprime  = xNewZeroInit<IssmDouble>(dim*numdof);
+                IssmDouble* Dlambda  = xNewZeroInit<IssmDouble>(dim*dim);
+                IssmDouble* C        = xNewZeroInit<IssmDouble>(dim*lnumdof);
+                IssmDouble* Cprime   = xNewZeroInit<IssmDouble>(dim*numdof);
+                IssmDouble* CtCUzawa = xNewZeroInit<IssmDouble>(numdof*lnumdof);
 
                 delete gauss;
-                Gauss* gauss=element->NewGaussBase(5);
+                gauss = element->NewGaussBase(5);
                 for(int ig=gauss->begin();ig<gauss->end();ig++){
                         gauss->GaussPoint(ig);
 
                         element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
                         this->GetCFS(C,element,dim,xyz_list,gauss);
                         this->GetCFSprime(Cprime,element,dim,xyz_list,gauss);
-                        for(i=0;i<dim;i++)   Dlambda[i*epssize+i] = gauss->weight*Jdet*sqrt(normal[i]*normal[i])*sqrt(rl);
+                        for(i=0;i<dim;i++) Dlambda[i*dim+i] = gauss->weight*Jdet*sqrt(normal[i]*normal[i])*sqrt(rl);
                         TripleMultiply(C,dim,lnumdof,1,
                                                 Dlambda,dim,dim,0,
                                                 Cprime,dim,numdof,0,
                                                 CtCUzawa,1);
                 }
+
+                /*The sigma naugmentation should not be transformed*/
+                MatrixMultiply(CtCUzawa,lnumdof,numdof,1,
+                                        CtCUzawa,lnumdof,numdof,0,
+                                        &Ke->values[0],1);
+
                 /*Delete base part*/
-                xDelete<IssmDouble>(xyz_list_base);
                 xDelete<IssmDouble>(Dlambda);
                 xDelete<IssmDouble>(C);
                 xDelete<IssmDouble>(Cprime);
+                xDelete<IssmDouble>(CtCUzawa);
+                xDelete<IssmDouble>(xyz_list_base);
         }
 
         /*Transform Coordinate System*/
@@ -2998 +2997 @@
                                 BtBUzawa,pnumdof,numdof,0,
                                 &Ke->values[0],1);
 
-        /*The sigma naugmentation should not be transformed*/
-        MatrixMultiply(CtCUzawa,lnumdof,numdof,1,
-                                CtCUzawa,lnumdof,numdof,0,
-                                &Ke->values[0],1);
-
         /*Clean up and return*/
         delete gauss;
         xDelete<IssmDouble>(xyz_list);
@@ -3012 +3006 @@
         xDelete<IssmDouble>(BprimeU);
         xDelete<IssmDouble>(BU);
         xDelete<IssmDouble>(BtBUzawa);
-        xDelete<IssmDouble>(Cprime);
-        xDelete<IssmDouble>(C);
-        xDelete<IssmDouble>(CtCUzawa);
         xDelete<int>(cs_list);
         return Ke;
 }/*}}}*/
@@ -3297 +3288 @@
 }/*}}}*/
 ElementVector* StressbalanceAnalysis::CreatePVectorFS(Element* element){/*{{{*/
 
+        ElementVector* pe = NULL;
+
+        ElementVector* pe1=CreatePVectorFSViscous(element);
+        ElementVector* pe2=CreatePVectorFSFriction(element);
+        ElementVector* pe3=CreatePVectorFSStress(element);
+        pe =new ElementVector(pe1,pe2,pe3);
+        delete pe1;
+        delete pe2;
+        delete pe3;
+
+        /*clean-up and return*/
+        return pe;
+}/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorFSViscous(Element* element){/*{{{*/
+
         int         i,dim,fe_FS;
         IssmDouble  x_coord,y_coord,z_coord;
         IssmDouble  Jdet,forcex,forcey,forcez;
@@ -3370 +3376 @@
         }
         return pe;
 }/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorFSFriction(Element* element){/*{{{*/
+
+        if(!element->IsOnBase()) return NULL;
+
+        /*Intermediaries*/
+        int         dim;
+        IssmDouble  alpha2,Jdet;
+        IssmDouble  bed_normal[3];
+        IssmDouble *xyz_list_base = NULL;
+        Gauss*      gauss         = NULL;
+
+        /*Get problem dimension*/
+        element->FindParam(&dim,DomainDimensionEnum);
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int vnumnodes = element->NumberofNodesVelocity();
+
+        /*Initialize Element matrix and vectors*/
+        ElementVector* pe = element->NewElementVector(FSvelocityEnum);
+        IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
+
+        /*Retrieve all inputs and parameters*/
+        element->GetVerticesCoordinatesBase(&xyz_list_base);
+        Input*  alpha2_input=element->GetInput(FrictionCoefficientEnum); _assert_(alpha2_input);
+
+        /* Start  looping on the number of gaussian points: */
+        gauss=element->NewGaussBase(3);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+                gauss->GaussPoint(ig);
+
+                alpha2_input->GetInputValue(&alpha2, gauss);
+                element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
+                element->NodalFunctionsVelocity(vbasis,gauss);
+                element->NormalBase(&bed_normal[0],xyz_list_base);
+
+                for(int i=0;i<vnumnodes;i++){
+                        pe->values[i*dim+0] += - alpha2*gauss->weight*Jdet*vbasis[i]*bed_normal[1];
+                        pe->values[i*dim+1] += alpha2*gauss->weight*Jdet*vbasis[i]*bed_normal[0];
+                        if(dim==3){
+                                pe->values[i*dim+2]+= alpha2*gauss->weight*Jdet*vbasis[i];
+                        }
+                }
+
+        }
+
+        /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/
+
+        /*Clean up and return*/
+        delete gauss;
+        xDelete<IssmDouble>(xyz_list_base);
+        xDelete<IssmDouble>(vbasis);
+        return pe;
+}/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorFSStress(Element* element){/*{{{*/
+
+        if(!element->IsOnBase()) return NULL;
+
+        /*Intermediaries*/
+        int         dim;
+        IssmDouble  sigmann,sigmant,Jdet,bedslope,beta;
+        IssmDouble *xyz_list_base = NULL;
+        Gauss*      gauss         = NULL;
+
+        /*Get problem dimension*/
+        element->FindParam(&dim,DomainDimensionEnum);
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int vnumnodes = element->NumberofNodesVelocity();
+
+        /*Initialize Element matrix and vectors*/
+        ElementVector* pe = element->NewElementVector(FSvelocityEnum);
+        IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
+
+        /*Retrieve all inputs and parameters*/
+        element->GetVerticesCoordinatesBase(&xyz_list_base);
+        Input*  sigmann_input=element->GetInput(VzEnum); _assert_(sigmann_input);
+        Input*  sigmant_input=element->GetInput(TemperatureEnum); _assert_(sigmant_input);
+        Input*  bedslope_input=element->GetInput(BedSlopeXEnum);     _assert_(bedslope_input);
+
+        /* Start  looping on the number of gaussian points: */
+        gauss=element->NewGaussBase(3);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+                gauss->GaussPoint(ig);
+
+                sigmann_input->GetInputValue(&sigmann, gauss);
+                sigmant_input->GetInputValue(&sigmant, gauss);
+                bedslope_input->GetInputValue(&bedslope, gauss);
+                element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
+                element->NodalFunctionsVelocity(vbasis,gauss);
+
+                beta=sqrt(1+bedslope*bedslope);
+                for(int i=0;i<vnumnodes;i++){
+                        pe->values[i*dim+0] += - (1./beta)*(-bedslope*sigmann + sigmant)*gauss->weight*Jdet*vbasis[i];
+                        pe->values[i*dim+1] += - (1./beta)*(sigmann + bedslope*sigmant)*gauss->weight*Jdet*vbasis[i];
+                        if(dim==3){
+                                //pe->values[i*dim+2]+= alpha2*gauss->weight*Jdet*vbasis[i];
+                                _error_("3d not supported yet");
+                        }
+                }
+
+        }
+
+        /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/
+
+        /*Clean up and return*/
+        delete gauss;
+        xDelete<IssmDouble>(xyz_list_base);
+        xDelete<IssmDouble>(vbasis);
+        return pe;
+}/*}}}*/
 #else
 ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSFriction(Element* element){/*{{{*/
 
@@ -3501 +3617 @@
         /*clean-up and return*/
         return pe;
 }/*}}}*/
-#endif
 ElementVector* StressbalanceAnalysis::CreatePVectorFSViscous(Element* element){/*{{{*/
 
         int         i,dim;
@@ -3571 +3686 @@
         xDelete<IssmDouble>(xyz_list);
         return pe;
 }/*}}}*/
+#endif
 ElementVector* StressbalanceAnalysis::CreatePVectorFSViscousXTH(Element* element){/*{{{*/
 
         int         i,tausize,dim;
@@ -3771 +3887 @@
         element->FindParam(&r,AugmentedLagrangianREnum);
         element->GetVerticesCoordinates(&xyz_list);
 
-        /*Get d and tau*/
+        /*Get pressure and sigmann*/
         Input* pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input);
         Input* sigmann_input =element->GetInput(SigmaNNEnum);  _assert_(sigmann_input);
 
@@ -3786 +3902 @@
                 for(i=0;i<numnodes;i++){
                         pe->values[i*dim+0] += pressure*gauss->weight*Jdet*dbasis[0*numnodes+i];
                         pe->values[i*dim+1] += pressure*gauss->weight*Jdet*dbasis[1*numnodes+i];
-                        if(dim==3){
-                                pe->values[i*dim+2]+= pressure*gauss->weight*Jdet*dbasis[2*numnodes+i];
-                        }
+                        if(dim==3) pe->values[i*dim+2]+= pressure*gauss->weight*Jdet*dbasis[2*numnodes+i];
                 }
         }
 
         if(element->IsOnBase()){ 
-
                 IssmDouble   sigmann;
                 IssmDouble*  vbasis = xNew<IssmDouble>(numnodes);
 
@@ -3801 +3914 @@
                 element->NormalBase(&bed_normal[0],xyz_list_base);
 
                 delete gauss;
-                Gauss* gauss=element->NewGaussBase(5);
+                gauss=element->NewGaussBase(5);
                 for(int ig=gauss->begin();ig<gauss->end();ig++){
                         gauss->GaussPoint(ig);
 
@@ -3810 +3923 @@
                         sigmann_input->GetInputValue(&sigmann, gauss);
 
                         for(i=0;i<numnodes;i++){
-                                pe->values[i*dim+0] += - sigmann*bed_normal[0]*gauss->weight*Jdet*vbasis[i];
-                                pe->values[i*dim+1] += - sigmann*bed_normal[1]*gauss->weight*Jdet*vbasis[i];
-                                if(dim==3){
-                                        pe->values[i*dim+2]+= - sigmann*bed_normal[2]*gauss->weight*Jdet*vbasis[i];
-                                }
+                                pe->values[i*dim+0] += + sigmann*bed_normal[0]*gauss->weight*Jdet*vbasis[i];
+                                pe->values[i*dim+1] += + sigmann*bed_normal[1]*gauss->weight*Jdet*vbasis[i];
+                                if(dim==3) pe->values[i*dim+2] += + sigmann*bed_normal[2]*gauss->weight*Jdet*vbasis[i];
                         }
                 }
                 xDelete<IssmDouble>(xyz_list_base);
@@ -3822 +3933 @@
         }
 
         /*Transform coordinate system*/
-        //element->TransformLoadVectorCoord(pe,cs_list); Do not transform pressure augmentation
+        //element->TransformLoadVectorCoord(pe,cs_list); Do not transform augmentation
 
         /*Clean up and return*/
         delete gauss;
@@ -4470 +4581 @@
          */
 
         /*Fetch number of nodes for this finite element*/
-        int lnumnodes;
-        if(dim==2) lnumnodes=3;
-        else lnumnodes=6;
-        //int pnumnodes = element->NumberofNodes(P1Enum);
+        int lnumnodes = element->GetNumberOfNodes(P2Enum);
 
         /*Get nodal functions derivatives*/
         IssmDouble* basis =xNew<IssmDouble>(lnumnodes);
-        element->NodalFunctions(basis,gauss);
+        element->NodalFunctionsP2(basis,gauss);
 
         /*Build B: */
         for(int i=0;i<lnumnodes;i++){
-                C[i*lnumnodes+0] = basis[i];
-                C[i*lnumnodes+1] = basis[i];
+                C[lnumnodes*0+i] = basis[i];
+                C[lnumnodes*1+i] = basis[i];
         }
 
         /*Clean up*/
@@ -4501 +4609 @@
          */
 
         /*Fetch number of nodes for this finite element*/
-        int vnumnodes = element->GetNumberOfNodes();
+        int vnumnodes = element->NumberofNodesVelocity();
         int vnumdof   = vnumnodes*dim;
 
         IssmDouble* vbasis=xNew<IssmDouble>(vnumnodes);