Changeset 22523

Timestamp:

03/11/18 14:38:34 (7 years ago)

Author:

Mathieu Morlighem

Message:

CHG: removing TripleMultiply where needed

Location:

issm/trunk-jpl/src/c/analyses

Files:

• LevelsetAnalysis.cpp (modified) (7 diffs)
• ThermalAnalysis.cpp (modified) (1 diff)

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

@@ -156 +156 @@
         /*Intermediaries */
         int  stabilization,dim, domaintype, calvinglaw;
-        int i, row, col;
+        int i,j,k,row, col;
         IssmDouble kappa;
         IssmDouble Jdet, dt, D_scalar;
         IssmDouble h,hx,hy,hz;
-        IssmDouble vel;
+        IssmDouble vel,v[3],w[3],c[3],m[3],dlsf[3];
         IssmDouble norm_dlsf, norm_calving, calvingrate, meltingrate, groundedice;
         IssmDouble calvingmax, calvinghaf, heaviside, haf_eps;
@@ -184 +184 @@
         ElementMatrix* Ke       = basalelement->NewElementMatrix();
         IssmDouble*    basis    = xNew<IssmDouble>(numnodes);
-        IssmDouble*    B        = xNew<IssmDouble>(dim*numnodes);
-        IssmDouble*    Bprime   = xNew<IssmDouble>(dim*numnodes);
-        IssmDouble*    D        = xNew<IssmDouble>(dim*dim);
-        IssmDouble*    v        = xNew<IssmDouble>(dim);
-        IssmDouble*    w        = xNew<IssmDouble>(dim);
-        IssmDouble*    c        = xNewZeroInit<IssmDouble>(dim);
-        IssmDouble*    m        = xNewZeroInit<IssmDouble>(dim);
-        IssmDouble*    dlsf     = xNew<IssmDouble>(dim);
+        IssmDouble*    dbasis   = xNew<IssmDouble>(2*numnodes);
+        IssmDouble*    Bprime = NULL;
+        if(stabilization==2){
+                Bprime   = xNew<IssmDouble>(dim*numnodes);
+        }
 
         /*Retrieve all inputs and parameters*/
@@ -283 +280 @@
 
                 basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                basalelement->NodalFunctions(basis,gauss);
+                basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
                 D_scalar=gauss->weight*Jdet;
 
                 /* Transient */
                 if(dt!=0.){
-                        basalelement->NodalFunctions(basis,gauss);
-                        TripleMultiply(basis,numnodes,1,0,
-                                                &D_scalar,1,1,0,
-                                                basis,1,numnodes,0,
-                                                &Ke->values[0],1);
-                        D_scalar*=dt;
+                        for(i=0;i<numnodes;i++){
+                                for(j=0;j<numnodes;j++){
+                                        Ke->values[i*numnodes+j] += D_scalar*basis[j]*basis[i];
+                                }
+                        }
+                        D_scalar=D_scalar*dt;
                 }
 
                 /* Advection */
-                GetB(B,basalelement,xyz_list,gauss); 
-                GetBprime(Bprime,basalelement,xyz_list,gauss); 
                 vx_input->GetInputValue(&v[0],gauss);
                 vy_input->GetInputValue(&v[1],gauss); 
@@ -458 +455 @@
 
                 /*Compute D*/
-                for(row=0;row<dim;row++){
-                        for(col=0;col<dim;col++){
-                                if(row==col)
-                                 D[row*dim+col]=D_scalar*w[row];
-                                else
-                                 D[row*dim+col]=0.;
+                for(i=0;i<numnodes;i++){
+                        for(j=0;j<numnodes;j++){
+                                for(k=0;k<dim;k++){
+                                        Ke->values[i*numnodes+j] += D_scalar*w[k]*dbasis[k*numnodes+j]*basis[i];
+                                }
                         }
                 }
-
-                TripleMultiply(B,dim,numnodes,1,
-                                        D,dim,dim,0,
-                                        Bprime,dim,numnodes,0,
-                                        &Ke->values[0],1);
 
                 /* Stabilization */
@@ -485 +476 @@
                                 h=sqrt( pow(hx*w[0]/vel,2) + pow(hy*w[1]/vel,2) ); 
                                 kappa=h*vel/2.;
-                                for(row=0;row<dim;row++)
-                                        for(col=0;col<dim;col++)
-                                        if(row==col)
-                                                D[row*dim+col]=D_scalar*kappa;
-                                        else
-                                                D[row*dim+col]=0.;
-
-                                TripleMultiply(Bprime,dim,numnodes,1,
-                                                        D,dim,dim,0,
-                                                        Bprime,dim,numnodes,0,
-                                                        &Ke->values[0],1);
+                                for(i=0;i<numnodes;i++){
+                                        for(j=0;j<numnodes;j++){
+                                                for(k=0;k<dim;k++){
+                                                        Ke->values[i*numnodes+j] += D_scalar*kappa*dbasis[k*numnodes+j]*dbasis[k*numnodes+i];
+                                                }
+                                        }
+                                }
                                 break;  
                         case 2:
@@ -501 +488 @@
                                 basalelement->ElementSizes(&hx,&hy,&hz);
                                 h=sqrt( pow(hx*w[0]/vel,2) + pow(hy*w[1]/vel,2) );
+                                IssmDouble D[3*3];
                                 for(row=0;row<dim;row++) 
                                         for(col=0;col<dim;col++) 
                                                 D[row*dim+col] = D_scalar*h/(2.*vel)*w[row]*w[col];
+                                GetBprime(Bprime,basalelement,xyz_list,gauss);
 
                                 TripleMultiply(Bprime,dim,numnodes,1,
-                                                        D,dim,dim,0,
+                                                        &D[0],dim,dim,0,
                                                         Bprime,dim,numnodes,0,
                                                         &Ke->values[0],1);
@@ -518 +507 @@
         xDelete<IssmDouble>(xyz_list);
         xDelete<IssmDouble>(basis);
-        xDelete<IssmDouble>(B);
-        xDelete<IssmDouble>(D);
+        xDelete<IssmDouble>(dbasis);
         xDelete<IssmDouble>(Bprime);
-        xDelete<IssmDouble>(v);
-        xDelete<IssmDouble>(w);
-        xDelete<IssmDouble>(c);
-        xDelete<IssmDouble>(m);
-        xDelete<IssmDouble>(dlsf);
         delete gauss;
         if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};

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

@@ -461 +461 @@
                 }
                 else if(stabilization==2){
-                        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
                         tau_parameter=element->StabilizationParameter(u-um,v-vm,w-wm,diameter,kappa);
                         for(int i=0;i<numnodes;i++){