Changeset 18281

Timestamp:

07/22/14 12:50:33 (11 years ago)

Author:

seroussi

Message:

NEW: added second augmentation

File:

• issm/trunk-jpl/src/c/analyses/UzawaPressureAnalysis.cpp (modified) (3 diffs)

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

@@ -31 +31 @@
         if(iomodel->domaintype==Domain3DEnum) iomodel->FetchDataToInput(elements,VzEnum,0.);
         iomodel->FetchDataToInput(elements,PressureEnum,0.);
+        iomodel->FetchDataToInput(elements,SigmaNNEnum,0.);
 }/*}}}*/
 void UzawaPressureAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
@@ -171 +172 @@
         int        *doflist       = NULL;
         IssmDouble rholambda,un,vx,vy,vz,sigmann;
-        IssmDouble bed_normal[3];
         IssmDouble *xyz_list_base = NULL;
 
         /*Fetch number of nodes and dof for this finite element*/
-        int numnodes       = element->GetNumberOfNodes();
-        //int numnodessigma  = element->NumberofNodes(P2Enum);
-        int numnodessigma;
+        int numnodes      = element->GetNumberOfNodes();
+        int numnodessigma = element->GetNumberOfNodes(P2Enum);
 
         element->FindParam(&dim,DomainDimensionEnum);
-        if(dim==2) numnodessigma=3;
-        else numnodessigma=6;
 
         /*Fetch dof list and allocate solution vector*/
         element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-        IssmDouble* values    = 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);
+        IssmDouble* values        = xNew<IssmDouble>(numnodes);
+        IssmDouble* valueslambda  = xNewZeroInit<IssmDouble>(numnodessigma);
+        IssmDouble* pressure      = xNew<IssmDouble>(numnodes);
+        Input* sigmann_input      = element->GetInput(SigmaNNEnum); _assert_(sigmann_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);}
         element->GetInputListOnNodes(&pressure[0],PressureEnum);
 
         /*Update pressure enum first*/
         for(int i=0;i<numnodes;i++){
-                values[i]   = pressure[i] + solution[doflist[i]];
+                values[i]  = pressure[i] + solution[doflist[i]];
         }
         element->AddInput(PressureEnum,values,element->GetElementType());
@@ -202 +200 @@
         /*Now compute sigmann if on base*/
         if(element->IsOnBase()){ 
-
+                if(dim==3) _error_("not implemented yet");
+
+                int baselist[3];
+                int onbase=0;
+                IssmDouble  Jdet;
+                IssmDouble bed_normal[3];
+                IssmDouble  Jlambda[3][3]  = {0.0};
+                IssmDouble  Cuk[3]         = {0.0};
+                IssmDouble  deltalambda[3] = {0.0};
+                IssmDouble* vertexonbase  = xNew<IssmDouble>(numnodessigma);
+                Input* vertexonbase_input = element->GetInput(MeshVertexonbaseEnum); _assert_(vertexonbase_input);
+                Gauss* gauss = element->NewGauss();
+
+                IssmDouble* basis = xNewZeroInit<IssmDouble>(numnodessigma);
                 element->GetVerticesCoordinatesBase(&xyz_list_base);
                 element->NormalBase(&bed_normal[0],xyz_list_base);
                 element->FindParam(&rholambda,AugmentedLagrangianRholambdaEnum);
 
-                Gauss* gauss = element->NewGauss();
                 for(int i=0;i<numnodessigma;i++){
                         gauss->GaussNode(P2Enum,i);
-
-                        sigmann_input->GetInputValue(&sigmann, gauss);
+                        vertexonbase_input->GetInputValue(&vertexonbase[i], gauss);
+                        if(vertexonbase[i]==1){ 
+                                baselist[onbase]=i;
+                                onbase += 1;
+                        }
+                }
+                if(!onbase==3) _error_("basal nodes of element not found");
+
+                delete gauss;
+                gauss = element->NewGaussBase(3);
+                for(int ig=gauss->begin();ig<gauss->end();ig++){
+                        gauss->GaussPoint(ig);
+
+                        /*Compute Jlambda*/
+                        element->NodalFunctionsP2(basis,gauss);
+                        element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
+                        for(int i=0;i<3;i++){
+                                for(int j=0;j<3;j++){
+                                        Jlambda[i][j] += Jdet*gauss->weight*basis[baselist[i]]*basis[baselist[j]];
+                                }
+                        }
+
+                        /*Compute rho_lambd C u^k*/
                         vx_input->GetInputValue(&vx, gauss);
                         vy_input->GetInputValue(&vy, gauss);
                         un=bed_normal[0]*vx + bed_normal[1]*vy;
-                        if(dim==3){
-                           vz_input->GetInputValue(&vz, gauss);
-                                un = un + bed_normal[2]*vz;
-                        }
-                        values[i] = sigmann + rholambda*un;
-                }
+                        for(int i=0;i<3;i++) Cuk[i] += - un*rholambda*Jdet*gauss->weight*basis[baselist[i]];
+                }
+
+                /*Now update sigmann*/
+                Matrix3x3Solve(&deltalambda[0],&Jlambda[0][0],&Cuk[0]);
                 delete gauss;
+                gauss = element->NewGauss();
+                for(int i=0;i<3;i++){
+                        gauss->GaussNode(P2Enum,baselist[i]);
+                        sigmann_input->GetInputValue(&sigmann, gauss);
+                        valueslambda[baselist[i]] = sigmann + deltalambda[i];
+                }
+
+                delete gauss;
+                xDelete<IssmDouble>(vertexonbase);
                 xDelete<IssmDouble>(xyz_list_base);
-        }
-        element->AddInput(SigmaNNEnum,values,P2Enum);
+                xDelete<IssmDouble>(basis);
+        }
+        element->AddInput(SigmaNNEnum,valueslambda,P2Enum);
 
         /*Free ressources:*/
         xDelete<IssmDouble>(values);
+        xDelete<IssmDouble>(valueslambda);
         xDelete<IssmDouble>(pressure);
         xDelete<int>(doflist);
-
 }/*}}}*/
 void UzawaPressureAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/