Changeset 18235

Timestamp:

07/11/14 11:00:33 (11 years ago)

Author:

seroussi

Message:

NEW: trying to implement augmentation lambda

Location:

issm/trunk-jpl/src/c/analyses

Files:

• StressbalanceAnalysis.cpp (modified) (9 diffs)
• UzawaPressureAnalysis.cpp (modified) (2 diffs)

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

@@ -2881 +2881 @@
         int         i,dim,epssize;
         IssmDouble  r,rl,Jdet,viscosity,DU,DUl;
+        IssmDouble      normal[3];
         IssmDouble *xyz_list = NULL;
         IssmDouble *xyz_list_base = NULL;
@@ -2959 +2960 @@
                 element->FindParam(&rl,AugmentedLagrangianRlambdaEnum);
                 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);
 
                 delete gauss;
@@ -2966 +2973 @@
 
                         element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
-                        this->GetCFS(C,element,dim,xyz_list_base,gauss);
-                        this->GetCFSprime(Cprime,element,dim,xyz_list_base,gauss);
-                        DUl = gauss->weight*Jdet*sqrt(rl);
-                        TripleMultiply(C,1,lnumdof,1,
-                                                &DU,1,1,0,
-                                                Cprime,1,numdof,0,
+                        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);
+                        TripleMultiply(C,dim,lnumdof,1,
+                                                Dlambda,dim,dim,0,
+                                                Cprime,dim,numdof,0,
                                                 CtCUzawa,1);
                 }
+                /*Delete base part*/
+                xDelete<IssmDouble>(xyz_list_base);
+                xDelete<IssmDouble>(Dlambda);
+                xDelete<IssmDouble>(C);
+                xDelete<IssmDouble>(Cprime);
         }
 
@@ -3733 +3745 @@
         int         i,dim;
         IssmDouble  Jdet,r,pressure;
-        IssmDouble *xyz_list = NULL;
-        Gauss*      gauss    = NULL;
+        IssmDouble  bed_normal[3];
+        IssmDouble *xyz_list      = NULL;
+        IssmDouble *xyz_list_base = NULL;
+        Gauss*      gauss         = NULL;
 
         /*Get problem dimension*/
@@ -3757 +3771 @@
         /*Get d and tau*/
         Input* pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input);
+        Input* sigmann_input =element->GetInput(SigmaNNEnum);  _assert_(sigmann_input);
 
         gauss=element->NewGauss(5);
@@ -3773 +3788 @@
                         }
                 }
+        }
+
+        if(element->IsOnBase()){ 
+
+                IssmDouble   sigmann;
+                IssmDouble*  vbasis = xNew<IssmDouble>(numnodes);
+
+                element->GetVerticesCoordinatesBase(&xyz_list_base);
+                element->NormalBase(&bed_normal[0],xyz_list_base);
+
+                delete gauss;
+                Gauss* gauss=element->NewGaussBase(5);
+                for(int ig=gauss->begin();ig<gauss->end();ig++){
+                        gauss->GaussPoint(ig);
+
+                        element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
+                        element->NodalFunctionsVelocity(vbasis,gauss);
+                        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];
+                                }
+                        }
+                }
+                xDelete<IssmDouble>(xyz_list_base);
+                xDelete<IssmDouble>(vbasis);
         }
 
@@ -4410 +4454 @@
 }/*}}}*/
 void StressbalanceAnalysis::GetCFS(IssmDouble* C,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-        /*Compute C  matrix. C=[Cp1 Cp2 ...] where Cpi=phi_pi. 
+        /*Compute C  matrix. C=[Cp1 Cp2 ...] where:
+         *     Cpi=[phi phi]. 
          */
 
@@ -4421 +4466 @@
         /*Get nodal functions derivatives*/
         IssmDouble* basis =xNew<IssmDouble>(lnumnodes);
-        element->NodalFunctionsP1(basis,gauss);
+        element->NodalFunctions(basis,gauss);
 
         /*Build B: */
         for(int i=0;i<lnumnodes;i++){
-                C[i] = basis[i];
+                C[i*lnumnodes+0] = basis[i];
+                C[i*lnumnodes+1] = basis[i];
         }
 
@@ -4432 +4478 @@
 }/*}}}*/
 void StressbalanceAnalysis::GetCFSprime(IssmDouble* Cprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-        /*      Compute C'  matrix. C'=[C1' C2' C3'] 
-         *                      Ci' = [  dphi/dx   dphi/dy ]
+        /*      Compute C'  matrix. C'=[C1' C2' ...] 
+         *                      Ci' = [  phi  0  ]
+         *                            [   0  phi ]
          *
          *      In 3d
-         *         Ci=[ dh/dx   dh/dy    dh/dz  ]
+         *                      Ci' = [  phi  0   0  ]
+         *                            [   0  phi  0  ]
+         *                            [   0   0  phi ]
          *      where phi is the finiteelement function for node i.
          */
 
         /*Fetch number of nodes for this finite element*/
-        int lnumnodes = element->GetNumberOfNodes();
-
-        /*Get nodal functions derivatives*/
-        IssmDouble* dbasis=xNew<IssmDouble>(dim*lnumnodes);
-        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-
-        /*Build C_prime: */
-        if(dim==2){
-                for(int i=0;i<lnumnodes;i++){
-                        Cprime[dim*i+0] = dbasis[0*lnumnodes+i];
-                        Cprime[dim*i+1] = dbasis[1*lnumnodes+i];
+        int vnumnodes = element->GetNumberOfNodes();
+        int vnumdof   = vnumnodes*dim;
+
+        IssmDouble* vbasis=xNew<IssmDouble>(vnumnodes);
+        element->NodalFunctionsVelocity(vbasis,gauss);
+
+        /*Build B: */
+        if(dim==3){
+                for(int i=0;i<vnumnodes;i++){
+                        Cprime[vnumdof*0+3*i+0] = vbasis[i];
+                        Cprime[vnumdof*0+3*i+1] = 0.;
+                        Cprime[vnumdof*0+3*i+2] = 0.;
+
+                        Cprime[vnumdof*1+3*i+0] = 0.;
+                        Cprime[vnumdof*1+3*i+1] = vbasis[i];
+                        Cprime[vnumdof*1+3*i+2] = 0.;
+
+                        Cprime[vnumdof*2+3*i+0] = 0.;
+                        Cprime[vnumdof*2+3*i+1] = 0.;
+                        Cprime[vnumdof*2+3*i+2] = vbasis[i];
                 }
         }
         else{
-                for(int i=0;i<lnumnodes;i++){
-                        Cprime[dim*i+0] = dbasis[0*lnumnodes+i];
-                        Cprime[dim*i+1] = dbasis[1*lnumnodes+i];
-                        Cprime[dim*i+2] = dbasis[2*lnumnodes+i];
-                }
-        }
-
-        /*Clean up*/
-        xDelete<IssmDouble>(dbasis);
+                for(int i=0;i<vnumnodes;i++){
+                        Cprime[vnumdof*0+2*i+0] = vbasis[i];
+                        Cprime[vnumdof*0+2*i+1] = 0.;
+
+                        Cprime[vnumdof*1+2*i+0] = 0.;
+                        Cprime[vnumdof*1+2*i+1] = vbasis[i];
+                }
+        }
+
+        /*Clean-up*/
+        xDelete<IssmDouble>(vbasis);
 }/*}}}*/
 void StressbalanceAnalysis::GetSolutionFromInputsFS(Vector<IssmDouble>* solution,Element* element){/*{{{*/

issm/trunk-jpl/src/c/analyses/UzawaPressureAnalysis.cpp

@@ -12 +12 @@
 
         parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianRhopEnum));
+        parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianRholambdaEnum));
 }/*}}}*/
 void UzawaPressureAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
@@ -167 +168 @@
 void UzawaPressureAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
 
+        int        dim;
         int        *doflist   = NULL;
+        IssmDouble rholambda,un;
+        IssmDouble bed_normal[3];
 
         /*Fetch number of nodes and dof for this finite element*/
-        int numnodes = element->GetNumberOfNodes();
+        int numnodes       = element->GetNumberOfNodes();
+        int numnodessigma  = element->NumberofNodes(P2Enum);
+        element->FindParam(&dim,DomainDimensionEnum);
 
         /*Fetch dof list and allocate solution vector*/
         element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
         IssmDouble* values    = xNew<IssmDouble>(numnodes);
-        IssmDouble* pressure = xNew<IssmDouble>(numnodes);
+        IssmDouble* pressure  = xNew<IssmDouble>(numnodes);
+        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);}
+        Input* sigmann_input=element->GetInput(SigmaNNEnum); _assert_(vy_input);
         element->GetInputListOnNodes(&pressure[0],PressureEnum);
 
+        /*Update pressure enum first*/
         for(int i=0;i<numnodes;i++){
                 values[i]   = pressure[i] + solution[doflist[i]];
         }
-
         element->AddInput(PressureEnum,values,element->GetElementType());
 
+        /*Now compute sigmann if on base*/
+        element->GetInputListOnNodes(&sigmann[0],SigmaNNEnum);
+
+        if(element->IsOnBase()){ 
+
+                element->NormalBase(&bed_normal[0],xyz_list_tria);
+                element->FindParam(&rholambda,AugmentedLagrangianRholambdaEnum);
+
+                Gauss* gauss = element->NewGauss();
+                for(int i=0;i<numnodessigma;i++){
+                        gauss->GaussNode(P2Enum,i);
+
+                        vx_input->GetInputValue(&vx, gauss);
+                        vy_input->GetInputValue(&vy, gauss);
+                        un=bed_normal[0]*vx[i] + bed_normal[1]*vy[i];
+                        if(dim==3){
+                           vz_input->GetInputValue(&vz, gauss);
+                                un = un + bed_normal[2]*vz[i];
+                        }
+                        values[i] = sigmann[i] + rholambda*un;
+                }
+        }
+        element->AddInput(SigmaNNEnum,values,P2Enum));
+
         /*Free ressources:*/
+        delete gauss;
         xDelete<IssmDouble>(values);
+        xDelete<IssmDouble>(sigmann);
         xDelete<IssmDouble>(pressure);
         xDelete<int>(doflist);