Changeset 24431

Timestamp:

12/04/19 21:11:16 (5 years ago)

Author:

Mathieu Morlighem

Message:

CHG: trying to speed up the code

Location:

issm/trunk-jpl/src/c

Files:

• analyses/ExtrapolationAnalysis.cpp (modified) (2 diffs)
• analyses/ExtrapolationAnalysis.h (modified) (1 diff)
• analyses/MasstransportAnalysis.cpp (modified) (7 diffs)
• classes/Elements/Element.cpp (modified) (1 diff)
• classes/Elements/Element.h (modified) (1 diff)
• modules/SystemMatricesx/SystemMatricesx.cpp (modified) (2 diffs)

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

@@ -207 +207 @@
 ElementVector* ExtrapolationAnalysis::CreatePVector(Element* element){/*{{{*/
         return NULL;
-
-}/*}}}*/
-void           ExtrapolationAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss, int dim){/*{{{*/
-        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
-         * 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(NDOF1*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 d=0;d<dim;d++)
-                        B[numnodes*d+i] = basis[i];
-
-        /*Clean-up*/
-        xDelete<IssmDouble>(basis);
-}/*}}}*/
-void           ExtrapolationAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss, int dim){/*{{{*/
-        /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
-         * 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(NDOF2*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 d=0;d<dim;d++)
-                        Bprime[numnodes*d+i] = dbasis[d*numnodes+i];
-
-        /*Clean-up*/
-        xDelete<IssmDouble>(dbasis);
-
 }/*}}}*/
 void           ExtrapolationAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
@@ -288 +234 @@
         }
 }/*}}}*/
-int                             ExtrapolationAnalysis::GetExtrapolationCase(Element* element){/*{{{*/
+int            ExtrapolationAnalysis::GetExtrapolationCase(Element* element){/*{{{*/
 
         /* Get case of extrapolation, depending on domain quality, and extrapolation variable */

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

@@ -26 +26 @@
         ElementMatrix* CreateKMatrix(Element* element);
         ElementVector* CreatePVector(Element* element);
-        void           GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss, int dim);
-        void           GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss, int dim);
         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/MasstransportAnalysis.cpp

@@ -291 +291 @@
         IssmDouble xi,tau;
         IssmDouble dvx[2],dvy[2];
+        IssmDouble D[4];
         IssmDouble* xyz_list = NULL;
 
@@ -309 +310 @@
         IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
         IssmDouble*             dbasis = xNew<IssmDouble>(dim*numnodes);
-        IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
-        IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
-        IssmDouble*    D      = xNewZeroInit<IssmDouble>(dim*dim);
 
         /*Retrieve all inputs and parameters*/
@@ -333 +331 @@
                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
                 element->NodalFunctions(basis,gauss);
-
+                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+                /*Transient term*/
+                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[i]*basis[j];
+                        }
+                }
+
+                /*Advection terms*/
                 vxaverage_input->GetInputValue(&vx,gauss);
                 vxaverage_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+                D_scalar=dt*gauss->weight*Jdet;
                 if(dim==2){
                         vyaverage_input->GetInputValue(&vy,gauss);
                         vyaverage_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
-                }
-
-                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,element,dim,xyz_list,gauss);
-                GetBprime(Bprime,element,dim,xyz_list,gauss);
-
-                dvxdx=dvx[0];
-                if(dim==2) dvydy=dvy[1];
-                D_scalar=dt*gauss->weight*Jdet;
-
-                D[0*dim+0]=D_scalar*dvxdx;
-                if(dim==2) D[1*dim+1]=D_scalar*dvydy;
-
-                TripleMultiply(B,dim,numnodes,1,
-                                        D,dim,dim,0,
-                                        B,dim,numnodes,0,
-                                        &Ke->values[0],1);
-
-                D[0*dim+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);
-
+                        dvxdx=dvx[0];
+                        dvydy=dvy[1];
+                        for(int i=0;i<numnodes;i++){
+                                for(int j=0;j<numnodes;j++){
+                                        /*\phi_i \phi_j \nabla\cdot v*/
+                                        Ke->values[i*numnodes+j] += D_scalar*basis[i]*basis[j]*(dvxdx+dvydy);
+                                        /*\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]);
+                                }
+                        }
+                }
+                else{
+                        dvxdx=dvx[0];
+                        for(int i=0;i<numnodes;i++){
+                                for(int j=0;j<numnodes;j++){
+                                        Ke->values[i*numnodes+j] += D_scalar*dvxdx*dbasis[0*numnodes+i]*dbasis[0*numnodes+j];
+                                        Ke->values[i*numnodes+j] += D_scalar*vx*dbasis[0*numnodes+j]*basis[i];
+                                }
+                        }
+                }
+
+                for(int i=0;i<4;i++) D[i]=0.;
                 switch(stabilization){
                         case 0:
@@ -383 +383 @@
                         case 2:
                                 /*Streamline upwinding*/
-                                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
                                 vxaverage_input->GetInputAverage(&vx);
                                 if(dim==1){
@@ -397 +396 @@
                                 /*SUPG*/
                                 if(dim!=2) _error_("Stabilization "<<stabilization<<" not supported yet for dim != 2");
-                                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
                                 vxaverage_input->GetInputAverage(&vx);
                                 vyaverage_input->GetInputAverage(&vy);
@@ -418 +416 @@
                         }
 
-                        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] += D_scalar*h/(2.*vel)*dbasis[0*numnodes+i] *D[0]*dbasis[0*numnodes+j];
+                                        }
+                                }
+                        }
                 }
                 if(stabilization==2){
@@ -478 +489 @@
         xDelete<IssmDouble>(basis);
         xDelete<IssmDouble>(dbasis);
-        xDelete<IssmDouble>(B);
-        xDelete<IssmDouble>(Bprime);
-        xDelete<IssmDouble>(D);
         delete gauss;
         return Ke;

issm/trunk-jpl/src/c/classes/Elements/Element.cpp

@@ -47 +47 @@
 
 /*Other*/
+bool       Element::AnyFSet(){/*{{{*/
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int numnodes = this->GetNumberOfNodes();
+
+        for(int i=0;i<numnodes;i++){
+                if(nodes[i]->fsize) return true;
+        }
+        return false;
+}/*}}}*/
 void       Element::ComputeLambdaS(){/*{{{*/
 

issm/trunk-jpl/src/c/classes/Elements/Element.h

@@ -63 +63 @@
                 /*bool               AllActive(void);*/
                 /*bool               AnyActive(void);*/
+                bool               AnyFSet(void);
                 void               ComputeLambdaS(void);
                 void               ComputeNewDamage();

issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp

@@ -43 +43 @@
                 for (i=0;i<femmodel->elements->Size();i++){
                         element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+                        if(!element->AnyFSet()) continue;
                         ElementMatrix* Ke = analysis->CreateKMatrix(element);
                         ElementVector* pe = analysis->CreatePVector(element);
@@ -76 +77 @@
         for (i=0;i<femmodel->elements->Size();i++){
                 element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+                if(!element->AnyFSet()) continue;
                 ElementMatrix* Ke = analysis->CreateKMatrix(element);
                 ElementVector* pe = analysis->CreatePVector(element);