Changeset 26205

Timestamp:

04/22/21 10:46:39 (4 years ago)

Author:

tsantos

Message:

CHG: working on mlho: computing vertically averaged velocities

File:

• issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (modified) (7 diffs)

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

@@ -2732 +2732 @@
 ElementMatrix* StressbalanceAnalysis::CreateKMatrixMLHO(Element* element){/*{{{*/
         
-        _error_("Mono Layer Higher-Order called, not fully tested. If you are sure in using it, comment this line.");
+        _error_("Mono Layer Higher-Order called, not fully tested. If you are sure of using it, comment this line.");
         
         /* Check if ice in element */
@@ -2914 +2914 @@
 
         /* Start  looping on the number of gaussian points: */
-        Gauss* gauss      = element->NewGauss(5);
+        Gauss* gauss      = element->NewGauss(2);
         Gauss* gauss_base = basalelement->NewGauss();
         while(gauss->next()){
@@ -3192 +3192 @@
         switch(domaintype){
                 case Domain2DhorizontalEnum:
+                        // itapopo FIXME use basalelement or element only here
                         basalelement = element;
                         break;
@@ -3216 +3217 @@
         IssmDouble* vsx       = xNew<IssmDouble>(numnodes);
         IssmDouble* vsy       = xNew<IssmDouble>(numnodes);
+        IssmDouble* vx        = xNew<IssmDouble>(numnodes);
+        IssmDouble* vy        = xNew<IssmDouble>(numnodes);
         IssmDouble* vz        = xNew<IssmDouble>(numnodes);
         IssmDouble* vel       = xNew<IssmDouble>(numnodes);
+        IssmDouble* n                    = xNew<IssmDouble>(numnodes);
         
         /*Use the dof list to index into the solution vector: */
@@ -3252 +3256 @@
         element->AddBasalInput(VxShearEnum,vshx,element->GetElementType());
         element->AddBasalInput(VyShearEnum,vshy,element->GetElementType());
-
-        /*Get Vz and compute vel (base)*/
-        basalelement->GetInputListOnNodes(&vz[0],VzEnum,0.);
-        for(i=0;i<numnodes;i++) vel[i]=sqrt(vbx[i]*vbx[i] + vby[i]*vby[i] + vz[i]*vz[i]); 
-
+        
         /*Add vx and vy as inputs to the tria element (base velocities): */
         element->AddBasalInput(VxBaseEnum,vbx,element->GetElementType());
         element->AddBasalInput(VyBaseEnum,vby,element->GetElementType());
-        element->AddBasalInput(VelEnum,vel,element->GetElementType()); //itapopo check this
         
         /*Add vx and vy as inputs to the tria element (surface velocities): */
@@ -3266 +3265 @@
         element->AddBasalInput(VySurfaceEnum,vsy,element->GetElementType());
 
+        /*Compute the vertically averaged velocities on each node*/
+        basalelement->GetInputListOnNodes(&n[0],MaterialsRheologyNEnum,0.); 
+   for(i=0;i<numnodes;i++){ //numnodes of the 2D mesh in which the MLHO is written
+                vx[i]=vbx[i]+vshx[i]*(n[i]+1)/(n[i]+2);
+                vy[i]=vby[i]+vshy[i]*(n[i]+1)/(n[i]+2);
+        }
+                
+        /*Get Vz and compute vel (vertically averaged vel)*/
+        basalelement->GetInputListOnNodes(&vz[0],VzEnum,0.); 
+        for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]); 
+
+        /*Add vx and vy as inputs to the tria element (vertically averaged velocities): */
+        element->AddBasalInput(VxEnum,vx,element->GetElementType());
+        element->AddBasalInput(VyEnum,vy,element->GetElementType());
+        element->AddBasalInput(VelEnum,vel,element->GetElementType()); 
+
         /*Free ressources:*/
         xDelete<IssmDouble>(vel);
+        xDelete<IssmDouble>(vx);
+        xDelete<IssmDouble>(vy);
         xDelete<IssmDouble>(vz);
         xDelete<IssmDouble>(vby);
@@ -3277 +3294 @@
         xDelete<IssmDouble>(values);
         xDelete<IssmDouble>(xyz_list);
+        xDelete<IssmDouble>(n);
         xDelete<int>(doflist);
         if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;};