Changeset 25265

Timestamp:

07/10/20 21:26:15 (5 years ago)

Author:

Mathieu Morlighem

Message:

CHG: removing TripleMultiply

Location:

issm/trunk-jpl/src/c/analyses

Files:

• DepthAverageAnalysis.cpp (modified) (3 diffs)
• DepthAverageAnalysis.h (modified) (1 diff)
• ThermalAnalysis.cpp (modified) (1 diff)

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

@@ -62 +62 @@
         /*Initialize Element vector and other vectors*/
         ElementMatrix* Ke     = element->NewElementMatrix();
-        IssmDouble*    B      = xNew<IssmDouble>(numnodes);
+        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
+        IssmDouble*    dbasis = xNew<IssmDouble>(dim*numnodes);
 
         /*Retrieve all inputs and parameters*/
@@ -73 +74 @@
 
                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-                GetB(B,element,dim,xyz_list,gauss);
-                D=gauss->weight*Jdet*dt;
+                element->NodalFunctions(basis,gauss);
+                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
                 /*vertical diffusion*/
-                TripleMultiply(B,1,numnodes,1,
-                                        &D,1,1,0,
-                                        B,1,numnodes,0,
-                                        &Ke->values[0],1);
+                D=gauss->weight*Jdet*dt;
+                for(int i=0;i<numnodes;i++){
+                        for(int j=0;j<numnodes;j++){
+                                Ke->values[i*numnodes+j] += D*dbasis[(dim-1)*numnodes+i]*dbasis[(dim-1)*numnodes+j];
+                        }
+                }
 
-                /*Next value*/
+                /*Next value (Transient)*/
                 D=gauss->weight*Jdet;
-                element->NodalFunctions(B,gauss);
-                TripleMultiply(B,numnodes,1,0,
-                                        &D,1,1,0,
-                                        B,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[j]*basis[i];
         } 
 
         /*Clean up and return*/
         xDelete<IssmDouble>(xyz_list);
-        xDelete<IssmDouble>(B);
+        xDelete<IssmDouble>(basis);
+        xDelete<IssmDouble>(dbasis);
         delete gauss;
         return Ke;
@@ -137 +137 @@
         return pe;
 }/*}}}*/
-void           DepthAverageAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-        /*      Compute B  matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz];
-                where hi is the interpolation function for node i.*/
-
-        /*Fetch number of nodes for this finite element*/
-        int numnodes = element->GetNumberOfNodes();
-
-        /*Get nodal functions derivatives*/
-        IssmDouble* dbasis=xNew<IssmDouble>(dim*numnodes);
-        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-
-        /*Build B: */
-        for(int i=0;i<numnodes;i++){
-                B[i] = dbasis[(dim-1)*numnodes+i];
-        }
-
-        /*Clean-up*/
-        xDelete<IssmDouble>(dbasis);
-}
-/*}}}*/
 void           DepthAverageAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
            _error_("not implemented yet");

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

@@ -26 +26 @@
                 ElementMatrix* CreateKMatrix(Element* element);
                 ElementVector* CreatePVector(Element* element);
-                void           GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
                 void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);

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

@@ -446 +446 @@
                 if(dt!=0.){
                         D_scalar=gauss->weight*Jdet;
-                        for(int i=0;i<numnodes;i++){
-                                for(int j=0;j<numnodes;j++){
-                                        Ke->values[i*numnodes+j] += D_scalar*basis[j]*basis[i];
-                                }
-                        }
+                        for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D_scalar*basis[j]*basis[i];
                         D_scalar=D_scalar*dt;
                 }