Changeset 25261

Timestamp:

07/10/20 20:29:59 (5 years ago)

Author:

Mathieu Morlighem

Message:

CHG: removing TripleMultiply

Location:

issm/trunk-jpl/src/c/analyses

Files:

• MeltingAnalysis.cpp (modified) (1 diff)
• ThermalAnalysis.cpp (modified) (5 diffs)
• ThermalAnalysis.h (modified) (1 diff)

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

@@ -109 +109 @@
                 D=latentheat/heatcapacity*gauss->weight*Jdet;
 
-                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];
         }
 

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

@@ -362 +362 @@
                 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];
         }
 
@@ -395 +391 @@
         IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
         IssmDouble*    dbasis = xNew<IssmDouble>(3*numnodes);
-        IssmDouble*    Bprime = xNew<IssmDouble>(3*numnodes);
         IssmDouble     K[3][3];
 
@@ -468 +463 @@
                         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]) 
+                                                                );
+                                }
+                        }
                 }
                 else if(stabilization==2){
@@ -513 +512 @@
         /*Clean up and return*/
         xDelete<IssmDouble>(xyz_list);
-        xDelete<IssmDouble>(Bprime);
         xDelete<IssmDouble>(basis);
         xDelete<IssmDouble>(dbasis);
@@ -744 +742 @@
         return pe;
 
-}/*}}}*/
-void           ThermalAnalysis::GetBAdvec(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*1.
-         * For node i, Bi' can be expressed in the actual coordinate system
-         * by:
-         *       Bi_advec =[ h ]
-         *                 [ h ]
-         *                 [ h ]
-         * where h is the interpolation function for node i.
-         *
-         * We assume B has been allocated already, of size: 3x(1*NUMNODESP1)
-         */
-
-        /*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];
-                B[numnodes*2+i] = basis[i];
-        }
-
-        /*Clean-up*/
-        xDelete<IssmDouble>(basis);
-}/*}}}*/
-void           ThermalAnalysis::GetBAdvecprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*1.
-         * For node i, Bi' can be expressed in the actual coordinate system
-         * by:
-         *       Biprime_advec=[ dh/dx ]
-         *                     [ dh/dy ]
-         *                     [ dh/dz ]
-         * where h is the interpolation function for node i.
-         *
-         * We assume B has been allocated already, of size: 3x(1*numnodes)
-         */
-
-        /*Fetch number of nodes for this finite element*/
-        int numnodes = element->GetNumberOfNodes();
-
-        /*Get nodal functions derivatives*/
-        IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
-        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-
-        /*Build B: */
-        for(int i=0;i<numnodes;i++){
-                B[numnodes*0+i] = dbasis[0*numnodes+i];
-                B[numnodes*1+i] = dbasis[1*numnodes+i];
-                B[numnodes*2+i] = dbasis[2*numnodes+i];
-        }
-
-        /*Clean-up*/
-        xDelete<IssmDouble>(dbasis);
-}/*}}}*/
-void           ThermalAnalysis::GetBConduct(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*1.
-         * For node i, Bi' can be expressed in the actual coordinate system
-         * by:
-         *       Bi_conduct=[ dh/dx ]
-         *                  [ dh/dy ]
-         *                  [ dh/dz ]
-         * where h is the interpolation function for node i.
-         *
-         * We assume B has been allocated already, of size: 3x(1*numnodes)
-         */
-
-        /*Fetch number of nodes for this finite element*/
-        int numnodes = element->GetNumberOfNodes();
-
-        /*Get nodal functions derivatives*/
-        IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
-        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-
-        /*Build B: */
-        for(int i=0;i<numnodes;i++){
-                B[numnodes*0+i] = dbasis[0*numnodes+i];
-                B[numnodes*1+i] = dbasis[1*numnodes+i];
-                B[numnodes*2+i] = dbasis[2*numnodes+i];
-        }
-
-        /*Clean-up*/
-        xDelete<IssmDouble>(dbasis);
 }/*}}}*/
 void           ThermalAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/

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

@@ -31 +31 @@
                 ElementVector* CreatePVectorShelf(Element* element);
                 ElementVector* CreatePVectorVolume(Element* element);
-                void           GetBAdvec(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
-                void           GetBAdvecprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
-                void           GetBConduct(IssmDouble* B,Element* element,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);