Changeset 18184

Timestamp:

06/25/14 12:14:09 (11 years ago)

Author:

seroussi

Message:

BUG: debugging LATaylorHood

Location:

issm/trunk-jpl/src/c/analyses

Files:

• StressbalanceAnalysis.cpp (modified) (5 diffs)
• StressbalanceAnalysis.h (modified) (1 diff)

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

@@ -2855 +2855 @@
         if(fe_FS==XTaylorHoodEnum)
          Ke1=CreateKMatrixFSViscousXTH(element);
-        else if(fe_FS==LATaylorHoodEnum){
+        else if(fe_FS==LATaylorHoodEnum)
          Ke1=CreateKMatrixFSViscousLATH(element);
-        }
         else
          Ke1=CreateKMatrixFSViscous(element);
@@ -2902 +2901 @@
         IssmDouble*    B        = xNew<IssmDouble>(epssize*numdof);
         IssmDouble*    Bprime   = xNew<IssmDouble>(epssize*numdof);
-        IssmDouble*    BtBUzawa = xNew<IssmDouble>(numdof*pnumdof);
+        IssmDouble*    BtBUzawa = xNewZeroInit<IssmDouble>(numdof*pnumdof);
         IssmDouble*    BU       = xNew<IssmDouble>(pnumdof);
         IssmDouble*    BprimeU  = xNew<IssmDouble>(numdof);
@@ -2920 +2919 @@
 
                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-                this->GetBFS(B,element,dim,xyz_list,gauss);
-                this->GetBFSprime(Bprime,element,dim,xyz_list,gauss);
+                this->GetBFSvel(B,element,dim,xyz_list,gauss);
+                this->GetBFSprimevel(Bprime,element,dim,xyz_list,gauss);
 
                 element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
@@ -4121 +4120 @@
         xDelete<IssmDouble>(pbasis);
 }/*}}}*/
+void StressbalanceAnalysis::GetBFSvel(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+        /*Compute B  matrix. B=[Bv1 Bv2 ... Bp1 Bp2 ...] where Bvi is of size 3*NDOF3. 
+         * For node i, Bvi can be expressed in the actual coordinate system
+         * by:     Bvi=[ dphi/dx          0        ]
+         *                                       [   0           dphi/dy     ]
+         *                                       [ 1/2*dphi/dy    1/2*dphi/dx]
+         *
+         *
+         *      In 3d:
+         *         Bvi=[ dh/dx          0             0      ]
+         *                                      [   0           dh/dy           0      ]
+         *                                      [   0             0           dh/dz    ]
+         *                                      [ 1/2*dh/dy    1/2*dh/dx        0      ]
+         *                                      [ 1/2*dh/dz       0         1/2*dh/dx  ]
+         *                                      [   0          1/2*dh/dz    1/2*dh/dy  ]
+         *
+         *      where phi is the finiteelement function for node i.
+         *      Same thing for Bb except the last column that does not exist.
+         */
+
+        /*Fetch number of nodes for this finite element*/
+        int vnumnodes = element->NumberofNodesVelocity();
+
+        /*Get nodal functions derivatives*/
+        IssmDouble* vdbasis=xNew<IssmDouble>(dim*vnumnodes);
+        element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);
+
+        /*Build B: */
+        if(dim==2){
+                for(int i=0;i<vnumnodes;i++){
+                        B[(dim*vnumnodes)*0+dim*i+0] = vdbasis[0*vnumnodes+i];
+                        B[(dim*vnumnodes)*0+dim*i+1] = 0.;
+                        B[(dim*vnumnodes)*1+dim*i+0] = 0.;
+                        B[(dim*vnumnodes)*1+dim*i+1] = vdbasis[1*vnumnodes+i];
+                        B[(dim*vnumnodes)*2+dim*i+0] = .5*vdbasis[1*vnumnodes+i];
+                        B[(dim*vnumnodes)*2+dim*i+1] = .5*vdbasis[0*vnumnodes+i];
+                }
+        }
+        else{
+                for(int i=0;i<vnumnodes;i++){
+                        B[(dim*vnumnodes)*0+dim*i+0] = vdbasis[0*vnumnodes+i];
+                        B[(dim*vnumnodes)*0+dim*i+1] = 0.;
+                        B[(dim*vnumnodes)*0+dim*i+2] = 0.;
+                        B[(dim*vnumnodes)*1+dim*i+0] = 0.;
+                        B[(dim*vnumnodes)*1+dim*i+1] = vdbasis[1*vnumnodes+i];
+                        B[(dim*vnumnodes)*1+dim*i+2] = 0.;
+                        B[(dim*vnumnodes)*2+dim*i+0] = 0.;
+                        B[(dim*vnumnodes)*2+dim*i+1] = 0.;
+                        B[(dim*vnumnodes)*2+dim*i+2] = vdbasis[2*vnumnodes+i];
+                        B[(dim*vnumnodes)*3+dim*i+0] = .5*vdbasis[1*vnumnodes+i];
+                        B[(dim*vnumnodes)*3+dim*i+1] = .5*vdbasis[0*vnumnodes+i];
+                        B[(dim*vnumnodes)*3+dim*i+2] = 0.;
+                        B[(dim*vnumnodes)*4+dim*i+0] = .5*vdbasis[2*vnumnodes+i];
+                        B[(dim*vnumnodes)*4+dim*i+1] = 0.;
+                        B[(dim*vnumnodes)*4+dim*i+2] = .5*vdbasis[0*vnumnodes+i];
+                        B[(dim*vnumnodes)*5+dim*i+0] = 0.;
+                        B[(dim*vnumnodes)*5+dim*i+1] = .5*vdbasis[2*vnumnodes+i];
+                        B[(dim*vnumnodes)*5+dim*i+2] = .5*vdbasis[1*vnumnodes+i];
+                }
+        }
+
+        /*Clean up*/
+        xDelete<IssmDouble>(vdbasis);
+}/*}}}*/
+void StressbalanceAnalysis::GetBFSprimevel(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+        /*      Compute B'  matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2. 
+         *      For node i, Bi' can be expressed in the actual coordinate system
+         *      by: 
+         *                      Bvi' = [  dphi/dx     0     ]
+         *                                       [     0      dphi/dy ]
+         *                                       [  dphi/dy   dphi/dx ]
+         *
+         *      In 3d
+         *         Bvi=[ dh/dx     0        0    ]
+         *                                      [   0      dh/dy      0    ]
+         *                                      [   0        0      dh/dz  ]
+         *                                      [ dh/dy    dh/dx      0    ]
+         *                                      [ dh/dz      0      dh/dx  ]
+         *                                      [   0      dh/dz    dh/dy  ]
+         *      where phi is the finiteelement function for node i.
+         *      In 3d:
+         */
+
+        /*Fetch number of nodes for this finite element*/
+        int vnumnodes = element->NumberofNodesVelocity();
+
+        /*Get nodal functions derivatives*/
+        IssmDouble* vdbasis=xNew<IssmDouble>(dim*vnumnodes);
+        element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);
+
+        /*Build B_prime: */
+        if(dim==2){
+                for(int i=0;i<vnumnodes;i++){
+                        Bprime[(dim*vnumnodes)*0+dim*i+0] = vdbasis[0*vnumnodes+i];
+                        Bprime[(dim*vnumnodes)*0+dim*i+1] = 0.;
+                        Bprime[(dim*vnumnodes)*1+dim*i+0] = 0.;
+                        Bprime[(dim*vnumnodes)*1+dim*i+1] = vdbasis[1*vnumnodes+i];
+                        Bprime[(dim*vnumnodes)*2+dim*i+0] = vdbasis[1*vnumnodes+i];
+                        Bprime[(dim*vnumnodes)*2+dim*i+1] = vdbasis[0*vnumnodes+i];
+                }
+        }
+        else{
+                for(int i=0;i<vnumnodes;i++){
+                        Bprime[(dim*vnumnodes)*0+dim*i+0] = vdbasis[0*vnumnodes+i];
+                        Bprime[(dim*vnumnodes)*0+dim*i+1] = 0.;
+                        Bprime[(dim*vnumnodes)*0+dim*i+2] = 0.;
+                        Bprime[(dim*vnumnodes)*1+dim*i+0] = 0.;
+                        Bprime[(dim*vnumnodes)*1+dim*i+1] = vdbasis[1*vnumnodes+i];
+                        Bprime[(dim*vnumnodes)*1+dim*i+2] = 0.;
+                        Bprime[(dim*vnumnodes)*2+dim*i+0] = 0.;
+                        Bprime[(dim*vnumnodes)*2+dim*i+1] = 0.;
+                        Bprime[(dim*vnumnodes)*2+dim*i+2] = vdbasis[2*vnumnodes+i];
+                        Bprime[(dim*vnumnodes)*3+dim*i+0] = vdbasis[1*vnumnodes+i];
+                        Bprime[(dim*vnumnodes)*3+dim*i+1] = vdbasis[0*vnumnodes+i];
+                        Bprime[(dim*vnumnodes)*3+dim*i+2] = 0.;
+                        Bprime[(dim*vnumnodes)*4+dim*i+0] = vdbasis[2*vnumnodes+i];
+                        Bprime[(dim*vnumnodes)*4+dim*i+1] = 0.;
+                        Bprime[(dim*vnumnodes)*4+dim*i+2] = vdbasis[0*vnumnodes+i];
+                        Bprime[(dim*vnumnodes)*5+dim*i+0] = 0.;
+                        Bprime[(dim*vnumnodes)*5+dim*i+1] = vdbasis[2*vnumnodes+i];
+                        Bprime[(dim*vnumnodes)*5+dim*i+2] = vdbasis[1*vnumnodes+i];
+                }
+        }
+
+        /*Clean up*/
+        xDelete<IssmDouble>(vdbasis);
+}/*}}}*/
 void StressbalanceAnalysis::GetBFSUzawa(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
         /*Compute B  matrix. B=[Bp1 Bp2 ...] where Bpi=phi_pi. 
@@ -4284 +4410 @@
         /*Add value to global vector*/
         solution->SetValues(vnumdof,vdoflist,vvalues,INS_VAL);
-        solution->SetValues(pnumdof,pdoflist,pvalues,INS_VAL);
+        if(pnumdof>0) solution->SetValues(pnumdof,pdoflist,pvalues,INS_VAL);
 
         /*Free ressources:*/

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

@@ -81 +81 @@
                 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);
+                void GetBFSprimevel(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
                 void GetBFSFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
                 void GetBFSUzawa(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);