source: issm/oecreview/Archive/24684-25833/ISSM-24917-24918.diff

../trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp

@@ -590 +590 @@
         ElementMatrix* Ke       = element->NewElementMatrix();
         IssmDouble*    basis    = xNew<IssmDouble>(numnodes);
         IssmDouble*    dbasis   = xNew<IssmDouble>(3*numnodes);
-        IssmDouble*    Bprime   = xNew<IssmDouble>(3*numnodes);
         IssmDouble     K[3][3];
 
         /*Retrieve all inputs and parameters*/
@@ -666 +665 @@
                         K[2][0]=h/(2.*vel)*fabs(vz*vx);  K[2][1]=h/(2.*vel)*fabs(vz*vy); K[2][2]=h/(2.*vel)*fabs(vz*vz);
                         for(int i=0;i<3;i++) for(int j=0;j<3;j++) K[i][j] = D_scalar*K[i][j];
 
-                        GetBAdvecprime(Bprime,element,xyz_list,gauss);
-                        TripleMultiply(Bprime,3,numnodes,1,
-                                                &K[0][0],3,3,0,
-                                                Bprime,3,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] *(K[0][0]*dbasis[0*numnodes+j] + K[0][1]*dbasis[1*numnodes+j]+ K[0][2]*dbasis[2*numnodes+j]) +
+                                                                dbasis[1*numnodes+i] *(K[1][0]*dbasis[0*numnodes+j] + K[1][1]*dbasis[1*numnodes+j]+ K[1][2]*dbasis[2*numnodes+j]) +
+                                                                dbasis[2*numnodes+i] *(K[2][0]*dbasis[0*numnodes+j] + K[2][1]*dbasis[1*numnodes+j]+ K[2][2]*dbasis[2*numnodes+j]) 
+                                                                );
+                                }
+                        }
                 }
                 /*SUPG*/
                 else if(stabilization==2){
@@ -714 +717 @@
         xDelete<IssmDouble>(xyz_list);
         xDelete<IssmDouble>(basis);
         xDelete<IssmDouble>(dbasis);
-        xDelete<IssmDouble>(Bprime);
         delete gauss;
         return Ke;
 }/*}}}*/
@@ -757 +759 @@
 
                 D=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_vel/(heatcapacity*rho_ice);
                 if(reCast<bool,IssmDouble>(dt)) D=dt*D;
-                TripleMultiply(basis,numnodes,1,0,
-                                        &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];
         }
 
         /*Clean up and return*/