Changeset 26531

Timestamp:

11/03/21 06:48:00 (3 years ago)

Author:

Cheng Gong

Message:

NEW: use P1xP4 for 3D MLHO

Location:

issm/trunk-jpl

Files:

• src/c/classes/Elements/Penta.cpp (modified) (2 diffs)
• src/c/classes/Elements/Tria.cpp (modified) (1 diff)
• test/Archives/Archive335.arch (modified) ( previous)

issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

@@ -179 +179 @@
                 else _error_("not implemented yet");
         }
-
-}
-/*}}}*/
-void       Penta::Recover3DMLHOInput(int targetVel_enum, int numnodes, IssmDouble* vb,  IssmDouble* vsh, IssmDouble* n, IssmDouble* H, IssmDouble* s){/*{{{*/
-   /* Recover the velocity acording to v=vb+(1-\zeta^{n+1})vsh, where \zeta=(s-z)/H
-    * The variables vb, vsh, n, H and s are all from the 2D horizontal mesh(Tria), with "numnodes" DOFs
-    * To project to penta the DOFs are doubled in size
-    *
-    */
-   _assert_(this->inputs);
-   if(!IsOnBase()) return;
-   else{
-      if(targetVel_enum==VxEnum || targetVel_enum==VyEnum){
-         IssmDouble vel[NUMVERTICES];
-         IssmDouble* xyz_list = NULL;
-         Penta* penta = this;
-         _assert_(NUMVERTICES-2*numnodes==0);
-
-         for(;;){
-            penta->GetVerticesCoordinates(&xyz_list);
-            for(int i=0;i<numnodes;i++) {
-               vel[i] = vb[i] + vsh[i]*(1.0-pow((s[i]-xyz_list[i*3+2])/H[i], (n[i]+1)));
-               vel[i+NUMVERTICES2D] = vb[i] + vsh[i]*(1-pow((s[i]-xyz_list[(i+NUMVERTICES2D)*3+2])/H[i], (n[i]+1)));
-            }
-                                xDelete<IssmDouble>(xyz_list);
-
-            /*Add to the bottom side of the element*/
-            penta->AddInput(targetVel_enum,&vel[0],P1Enum);
-            if (penta->IsOnSurface()) break;
-            penta=penta->GetUpperPenta(); _assert_(penta->Id()!=this->id);
-         }
-      }
-      else _error_("not implemented yet");
-   }
 
 }
@@ -3298 +3264 @@
 int        Penta::PressureInterpolation(void){/*{{{*/
         return PentaRef::PressureInterpolation(this->element_type);
+}
+/*}}}*/
+void       Penta::Recover3DMLHOInput(int targetVel_enum, int numnodes, IssmDouble* vb,  IssmDouble* vsh, IssmDouble* n, IssmDouble* H, IssmDouble* s){/*{{{*/
+   /* Recover the velocity acording to v=vb+(1-\zeta^{n+1})vsh, where \zeta=(s-z)/H
+    * The variables vb, vsh, n, H and s are all from the 2D horizontal mesh(Tria), with "numnodes" DOFs
+    * To project to penta the DOFs are doubled in size
+    *
+    */
+   _assert_(this->inputs);
+   if(!IsOnBase()) return;
+   else{
+      if(targetVel_enum==VxEnum || targetVel_enum==VyEnum){
+         IssmDouble vel[NUMVERTICES2D*5];
+         IssmDouble* xyz_list = NULL;
+                        IssmDouble zi;
+         Penta* penta = this;
+                        numnodes = penta->NumberofNodes(P1xP4Enum);
+         _assert_(NUMVERTICES2D*5-numnodes==0);
+                        GaussPenta gauss;
+
+         for(;;){
+            penta->GetVerticesCoordinates(&xyz_list);
+            for(int i=0;i<NUMVERTICES2D;i++) {
+                                        for (int j=0;j<5;j++){
+                                                /* Get z-coordinate of the current node */
+                                                gauss.GaussNode(P1xP4Enum, i+j*NUMVERTICES2D);
+                                                zi = this->GetZcoord(xyz_list, &gauss);
+                       vel[i+j*NUMVERTICES2D] = vb[i] + vsh[i]*(1.0-pow((s[i]-zi)/H[i], (n[i]+1)));
+                                        }
+            }
+                                xDelete<IssmDouble>(xyz_list);
+
+            /*Add to the bottom side of the element*/
+            penta->AddInput(targetVel_enum,&vel[0],P1xP4Enum);
+            if (penta->IsOnSurface()) break;
+            penta=penta->GetUpperPenta(); _assert_(penta->Id()!=this->id);
+         }
+      }
+      else _error_("not implemented yet");
+   }
+
 }
 /*}}}*/

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

@@ -2604 +2604 @@
                                 input->Serve(numindices,&indices[0]);
                                 break;
+                        case P1xP4Enum:
                         case P1DGEnum:
                         case P1bubbleEnum: