Changeset 25266 for issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp

Timestamp:

07/10/20 21:33:01 (5 years ago)

Author:

Mathieu Morlighem

Message:

CHG: removing TripleMultiply

File:

• issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp (modified) (6 diffs)

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

@@ -161 +161 @@
         ElementMatrix* Ke     = element->NewElementMatrix();
         IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
-        IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
-        IssmDouble*    Bprime = xNew<IssmDouble>(2*numnodes);
+        IssmDouble*             dbasis = xNew<IssmDouble>(2*numnodes);
         IssmDouble     D[2][2]={0.};
 
@@ -183 +182 @@
                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
                 element->NodalFunctions(basis,gauss);
+                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
                 vx_input->GetInputValue(&vx,gauss);
@@ -189 +189 @@
                 vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
 
+                /*Transient term*/
                 D_scalar=gauss->weight*Jdet;
-
-                TripleMultiply(basis,1,numnodes,1,
-                                        &D_scalar,1,1,0,
-                                        basis,1,numnodes,0,
-                                        Ke->values,1);
-
-                GetB(B,element,xyz_list,gauss);
-                GetBprime(Bprime,element,xyz_list,gauss);
-
+                for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D_scalar*basis[i]*basis[j];
+
+                /*Advection terms*/
                 dvxdx=dvx[0];
                 dvydy=dvy[1];
                 D_scalar=dt*gauss->weight*Jdet;
-
-                D[0][0]=D_scalar*dvxdx;
-                D[1][1]=D_scalar*dvydy;
-                TripleMultiply(B,2,numnodes,1,
-                                        &D[0][0],2,2,0,
-                                        B,2,numnodes,0,
-                                        &Ke->values[0],1);
-
-                D[0][0]=D_scalar*vx;
-                D[1][1]=D_scalar*vy;
-                TripleMultiply(B,2,numnodes,1,
-                                        &D[0][0],2,2,0,
-                                        Bprime,2,numnodes,0,
-                                        &Ke->values[0],1);
+                for(int i=0;i<numnodes;i++){
+                        for(int j=0;j<numnodes;j++){
+                                /*\phi_i \phi_j \nabla\cdot v*/
+                                Ke->values[i*numnodes+j] += D_scalar*basis[i]*basis[j]*(dvxdx+dvydy);
+                                /*\phi_i v\cdot\nabla\phi_j*/
+                                Ke->values[i*numnodes+j] += D_scalar*basis[i]*(vx*dbasis[0*numnodes+j] + vy*dbasis[1*numnodes+j]);
+                        }
+                }
 
                 /*Artificial diffusivity*/
@@ -223 +212 @@
                 D[0][1]=D_scalar*diffusivity*h/(2*vel)*vx*vy;
                 D[1][1]=D_scalar*diffusivity*h/(2*vel)*vy*vy;
-                TripleMultiply(Bprime,2,numnodes,1,
-                                        &D[0][0],2,2,0,
-                                        Bprime,2,numnodes,0,
-                                        &Ke->values[0],1);
+                for(int i=0;i<numnodes;i++){
+                        for(int j=0;j<numnodes;j++){
+                                Ke->values[i*numnodes+j] += (
+                                                        dbasis[0*numnodes+i] *(D[0][0]*dbasis[0*numnodes+j] + D[0][1]*dbasis[1*numnodes+j]) +
+                                                        dbasis[1*numnodes+i] *(D[1][0]*dbasis[0*numnodes+j] + D[1][1]*dbasis[1*numnodes+j]) 
+                                                        );
+                        }
+                }
         }
 
@@ -232 +225 @@
         xDelete<IssmDouble>(xyz_list);
         xDelete<IssmDouble>(basis);
-        xDelete<IssmDouble>(B);
-        xDelete<IssmDouble>(Bprime);
+        xDelete<IssmDouble>(dbasis);
         delete gauss;
         return Ke;
@@ -285 +277 @@
         delete gauss;
         return pe;
-}/*}}}*/
-void           HydrologyShreveAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*2. 
-         * For node i, Bi can be expressed in the actual coordinate system
-         * by: 
-         *       Bi=[ N ]
-         *          [ N ]
-         * where N is the finiteelement function for node i.
-         *
-         * We assume B_prog has been allocated already, of size: 2x(1*numnodes)
-         */
-
-        /*Fetch number of nodes for this finite element*/
-        int numnodes = element->GetNumberOfNodes();
-
-        /*Get nodal functions*/
-        IssmDouble* basis=xNew<IssmDouble>(numnodes);
-        element->NodalFunctions(basis,gauss);
-
-        /*Build B: */
-        for(int i=0;i<numnodes;i++){
-                B[numnodes*0+i] = basis[i];
-                B[numnodes*1+i] = basis[i];
-        }
-
-        /*Clean-up*/
-        xDelete<IssmDouble>(basis);
-}/*}}}*/
-void           HydrologyShreveAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-        /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*2. 
-         * For node i, Bi' can be expressed in the actual coordinate system
-         * by: 
-         *       Bi_prime=[ dN/dx ]
-         *                [ dN/dy ]
-         * where N is the finiteelement function for node i.
-         *
-         * We assume B' has been allocated already, of size: 3x(2*numnodes)
-         */
-
-        /*Fetch number of nodes for this finite element*/
-        int numnodes = element->GetNumberOfNodes();
-
-        /*Get nodal functions derivatives*/
-        IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
-        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-
-        /*Build B': */
-        for(int i=0;i<numnodes;i++){
-                Bprime[numnodes*0+i] = dbasis[0*numnodes+i];
-                Bprime[numnodes*1+i] = dbasis[1*numnodes+i];
-        }
-
-        /*Clean-up*/
-        xDelete<IssmDouble>(dbasis);
-
 }/*}}}*/
 void           HydrologyShreveAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/