Changeset 24936

Timestamp:

05/30/20 22:25:06 (5 years ago)

Author:

Mathieu Morlighem

Message:

CHG: removing TripleMultiply... inefficient

Location:

issm/trunk-jpl/src/c/analyses

Files:

• BalancethicknessAnalysis.cpp (modified) (5 diffs)
• StressbalanceVerticalAnalysis.cpp (modified) (2 diffs)

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

@@ -164 +164 @@
         /*Initialize Element vector and other vectors*/
         ElementMatrix* Ke     = element->NewElementMatrix();
-        IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
-        IssmDouble*    Bprime = xNew<IssmDouble>(2*numnodes);
+        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
+        IssmDouble*             dbasis = xNew<IssmDouble>(2*numnodes);
         IssmDouble     D[2][2];
 
@@ -190 +190 @@
 
                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-                GetB(B,element,xyz_list,gauss);
-                GetBprime(Bprime,element,xyz_list,gauss);
+                element->NodalFunctions(basis,gauss);
+                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
                 vxaverage_input->GetInputValue(&vx,gauss);
@@ -201 +201 @@
                 D_scalar=gauss->weight*Jdet;
 
-                D[0][0]=D_scalar*dvxdx;
-                D[0][1]=0.;
-                D[1][0]=0.;
-                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]);
+                        }
+                }
 
                 if(stabilization==1){
@@ -239 +232 @@
                         D[0][1]=D_scalar*D[0][1];
                         D[1][1]=D_scalar*D[1][1];
-                        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]) 
+                                                                );
+                                }
+                        }
                 }
         }
@@ -248 +246 @@
         /*Clean up and return*/
         xDelete<IssmDouble>(xyz_list);
-        xDelete<IssmDouble>(B);
-        xDelete<IssmDouble>(Bprime);
+        xDelete<IssmDouble>(basis);
+        xDelete<IssmDouble>(dbasis);
         delete gauss;
         return Ke;

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

@@ -216 +216 @@
                 D = -gauss->weight*Jdet*normal[2];
 
-                TripleMultiply( basis,1,numnodes,1,
-                                        &D,1,1,0,
-                                        basis,1,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] += D*basis[i]*basis[j];
         }
 
@@ -256 +253 @@
                 D = -gauss->weight*Jdet*normal[2];
 
-                TripleMultiply( basis,1,numnodes,1,
-                                        &D,1,1,0,
-                                        basis,1,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] += D*basis[i]*basis[j];
         }