Changeset 25262

Timestamp:

07/10/20 20:49:27 (5 years ago)

Author:

Mathieu Morlighem

Message:

CHG: removing TripleMultiply

Location:

issm/trunk-jpl/src/c/analyses

Files:

• BalancethicknessAnalysis.cpp (modified) (5 diffs)
• BalancethicknessAnalysis.h (modified) (1 diff)
• GLheightadvectionAnalysis.cpp (modified) (4 diffs)
• GLheightadvectionAnalysis.h (modified) (1 diff)

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

@@ -264 +264 @@
         /*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];
 
@@ -288 +288 @@
 
                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                element->NodalFunctions(basis,gauss);
+                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
                 vxaverage_input->GetInputValue(&vx,gauss);
                 vyaverage_input->GetInputValue(&vy,gauss);
@@ -294 +296 @@
                 D_scalar=gauss->weight*Jdet;
 
-                /*WARNING: B and Bprime are inverted compared to CG*/
-                GetB(Bprime,element,xyz_list,gauss);
-                GetBprime(B,element,xyz_list,gauss);
-
+                /*WARNING: inverted compared to CG*/
                 D_scalar = - gauss->weight*Jdet;
                 D[0][0]  = D_scalar*vx;
@@ -303 +302 @@
                 D[1][0]  = 0.;
                 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++){
+                                Ke->values[i*numnodes+j] += D_scalar*basis[i]*basis[j]*(dvxdx+dvydy);
+                                Ke->values[i*numnodes+j] += D_scalar*basis[j]*(vx*dbasis[0*numnodes+i] + vy*dbasis[1*numnodes+i]);
+                        }
+                }
         }
 
         /*Clean up and return*/
         xDelete<IssmDouble>(xyz_list);
-        xDelete<IssmDouble>(B);
-        xDelete<IssmDouble>(Bprime);
+        xDelete<IssmDouble>(basis);
+        xDelete<IssmDouble>(dbasis);
         delete gauss;
         return Ke;
@@ -420 +420 @@
         delete gauss;
         return pe;
-}/*}}}*/
-void           BalancethicknessAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*2. 
-         * For node i, Bi can be expressed in the actual coordinate system
-         * by: 
-         *       Bi=[ N ]
-         *          [ N ]
-         * where N is the finiteelement function for node i.
-         *
-         * We assume B_prog has been allocated already, of size: 2x(1*numnodes)
-         */
-
-        /*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];
-        }
-
-        /*Clean-up*/
-        xDelete<IssmDouble>(basis);
-}/*}}}*/
-void           BalancethicknessAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-        /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*2. 
-         * For node i, Bi' can be expressed in the actual coordinate system
-         * by: 
-         *       Bi_prime=[ dN/dx ]
-         *                [ dN/dy ]
-         * where N is the finiteelement function for node i.
-         *
-         * We assume B' has been allocated already, of size: 3x(2*numnodes)
-         */
-
-        /*Fetch number of nodes for this finite element*/
-        int numnodes = element->GetNumberOfNodes();
-
-        /*Get nodal functions derivatives*/
-        IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
-        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-
-        /*Build B': */
-        for(int i=0;i<numnodes;i++){
-                Bprime[numnodes*0+i] = dbasis[0*numnodes+i];
-                Bprime[numnodes*1+i] = dbasis[1*numnodes+i];
-        }
-
-        /*Clean-up*/
-        xDelete<IssmDouble>(dbasis);
-
 }/*}}}*/
 void           BalancethicknessAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/

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

@@ -30 +30 @@
                 ElementVector* CreatePVectorCG(Element* element);
                 ElementVector* CreatePVectorDG(Element* element);
-                void           GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
-                void           GetBprime(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);

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

@@ -109 +109 @@
                 basalelement->NodalFunctions(basis,gauss);
                 basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-                GetBprime(Bprime,basalelement,dim,xyz_list,gauss);
 
                 /*Get velocity*/
@@ -161 +160 @@
                         D[0][0]=D_scalar*h/(2.*vel)*fabs(vx*vx);  D[0][1]=D_scalar*h/(2.*vel)*fabs(vx*vy);
                         D[1][0]=D_scalar*h/(2.*vel)*fabs(vy*vx);  D[1][1]=D_scalar*h/(2.*vel)*fabs(vy*vy);
-                        TripleMultiply(Bprime,dim,numnodes,1,
-                                                &D[0][0],2,2,0,
-                                                Bprime,dim,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]) 
+                                                                );
+                                }
+                        }
                 }
         }
@@ -172 +175 @@
         xDelete<IssmDouble>(basis);
         xDelete<IssmDouble>(dbasis);
-        xDelete<IssmDouble>(Bprime);
         delete gauss;
         if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
@@ -235 +237 @@
         return;
 }/*}}}*/
-
-void           GLheightadvectionAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*2. 
-         * For node i, Bi can be expressed in the actual coordinate system
-         * by: 
-         *       Bi=[ N ]
-         *          [ N ]
-         * where N is the finiteelement function for node i.
-         *
-         * We assume B_prog has been allocated already, of size: 2x(1*numnodes)
-         */
-
-        /*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++){
-                for(int j=0;j<dim;j++){
-                        B[numnodes*j+i] = basis[i];
-                }
-        }
-
-        /*Clean-up*/
-        xDelete<IssmDouble>(basis);
-}/*}}}*/
-void           GLheightadvectionAnalysis::GetBprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-        /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*2. 
-         * For node i, Bi' can be expressed in the actual coordinate system
-         * by: 
-         *       Bi_prime=[ dN/dx ]
-         *                [ dN/dy ]
-         * where N is the finiteelement function for node i.
-         *
-         * We assume B' has been allocated already, of size: 3x(2*numnodes)
-         */
-
-        /*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++){
-                for(int j=0;j<dim;j++){
-                        Bprime[numnodes*j+i] = dbasis[j*numnodes+i];
-                }
-        }
-
-        /*Clean-up*/
-        xDelete<IssmDouble>(dbasis);
-
-}/*}}}*/

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

@@ -30 +30 @@
                 void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
                 void           UpdateConstraints(FemModel* femmodel);
-
-                void           GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
-                void           GetBprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
 };
 #endif