Changeset 25264 for issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp

Timestamp:

07/10/20 21:19:59 (5 years ago)

Author:

Mathieu Morlighem

Message:

CHG: removing TripleMultiply

File:

• issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp (modified) (5 diffs)

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

@@ -169 +169 @@
         ElementMatrix* Ke     = basalelement->NewElementMatrix(NoneApproximationEnum);
         IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
-        IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
-        IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
+        IssmDouble*    dbasis = xNew<IssmDouble>(dim*numnodes);
         IssmDouble*    D      = xNew<IssmDouble>(dim*dim);
 
@@ -189 +188 @@
                 basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
                 basalelement->NodalFunctions(basis,gauss);
+                basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
                 vx_input->GetInputValue(&vx,gauss);
                 if(dim==2) vy_input->GetInputValue(&vy,gauss);
 
+                /*Transient term*/
                 D_scalar=gauss->weight*Jdet;
-                TripleMultiply(basis,1,numnodes,1,
-                                        &D_scalar,1,1,0,
-                                        basis,1,numnodes,0,
-                                        &Ke->values[0],1);
-
-                GetB(B,basalelement,dim,xyz_list,gauss);
-                GetBprime(Bprime,basalelement,dim,xyz_list,gauss);
-
+                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];
+
+                /*Advection terms*/
                 D_scalar=dt*gauss->weight*Jdet;
                 for(int i=0;i<dim*dim;i++) D[i]=0.;
-                D[0] = D_scalar*vx;
-                if(dim==2) D[1*dim+1] = D_scalar*vy;
-
-                TripleMultiply(B,dim,numnodes,1,
-                                        D,dim,dim,0,
-                                        Bprime,dim,numnodes,0,
-                                        &Ke->values[0],1);
+
+                if(dim==1){
+                        /*\phi_i v\cdot\nabla\phi_j*/
+                        for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D_scalar*basis[i]*vx*dbasis[0*numnodes+j];
+                }
+                else{
+                        _assert_(dim==2);
+                        for(int i=0;i<numnodes;i++){
+                                for(int j=0;j<numnodes;j++){
+                                        /*\phi_i v\cdot\nabla\phi_j*/
+                                        Ke->values[i*numnodes+j] += D_scalar*basis[i]*(vx*dbasis[0*numnodes+j] + vy*dbasis[1*numnodes+j]);
+                                }
+                        }
+                }
 
                 if(stabilization==2){
@@ -241 +244 @@
                 if(stabilization==1 || stabilization==2){
                         for(int i=0;i<dim*dim;i++) D[i]=D_scalar*D[i];
-                        TripleMultiply(Bprime,dim,numnodes,1,
-                                                D,dim,dim,0,
-                                                Bprime,dim,numnodes,0,
-                                                &Ke->values[0],1);
+                        if(dim==2){
+                                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*dim+0]*dbasis[0*numnodes+j] + D[0*dim+1]*dbasis[1*numnodes+j]) +
+                                                                        dbasis[1*numnodes+i] *(D[1*dim+0]*dbasis[0*numnodes+j] + D[1*dim+1]*dbasis[1*numnodes+j]) 
+                                                                        );
+                                        }
+                                }
+                        }
+                        else{
+                                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]*dbasis[0*numnodes+j];
+                        }
                 }
         }
@@ -251 +263 @@
         xDelete<IssmDouble>(xyz_list);
         xDelete<IssmDouble>(basis);
-        xDelete<IssmDouble>(B);
-        xDelete<IssmDouble>(Bprime);
+        xDelete<IssmDouble>(dbasis);
         xDelete<IssmDouble>(D);
         delete gauss;
@@ -332 +343 @@
         if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
         return pe;
-
-}/*}}}*/
-void           FreeSurfaceBaseAnalysis::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           FreeSurfaceBaseAnalysis::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);
 
 }/*}}}*/