Changeset 22510

Timestamp:

03/09/18 11:32:45 (7 years ago)

Author:

Mathieu Morlighem

Message:

CHG: getting rid of TripleMultiply in HO (too much time lost with multiplications and xNew in GetB

File:

• issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (modified) (7 diffs)

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

@@ -2336 +2336 @@
 
         /*Initialize Element matrix and vectors*/
-        ElementMatrix* Ke = element->NewElementMatrix(HOApproximationEnum);
-        IssmDouble*    B  = xNew<IssmDouble>((dim-1)*numdof);
-        IssmDouble*    D  = xNewZeroInit<IssmDouble>((dim-1)*(dim-1));
+        ElementMatrix* Ke     = element->NewElementMatrix(HOApproximationEnum);
+        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
 
         /*Retrieve all inputs and parameters*/
@@ -2370 +2369 @@
                 }
 
-                this->GetBHOFriction(B,element,dim,xyz_list_base,gauss);
                 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
-                for(int i=0;i<dim-1;i++) D[i*(dim-1)+i]=alpha2*gauss->weight*Jdet;
-
-                TripleMultiply(B,dim-1,numdof,1,
-                                        D,dim-1,dim-1,0,
-                                        B,dim-1,numdof,0,
-                                        &Ke->values[0],1);
+                element->NodalFunctions(basis,gauss);
+
+                if(dim==3){
+                        for(int i=0;i<numnodes;i++){
+                                for(int j=0;j<numnodes;j++){
+                                        Ke->values[2*i*2*numnodes+2*j]       += alpha2*gauss->weight*Jdet*basis[i]*basis[j];
+                                        Ke->values[(2*i+1)*2*numnodes+2*j+1] += alpha2*gauss->weight*Jdet*basis[i]*basis[j];
+                                }
+                        }
+                }
+                else{
+                        for(int i=0;i<numnodes;i++){
+                                for(int j=0;j<numnodes;j++){
+                                        Ke->values[i*numnodes+j] += alpha2*gauss->weight*Jdet*basis[i]*basis[j];
+                                }
+                        }
+                }
         }
 
@@ -2387 +2396 @@
         delete friction;
         xDelete<IssmDouble>(xyz_list_base);
-        xDelete<IssmDouble>(B);
-        xDelete<IssmDouble>(D);
+        xDelete<IssmDouble>(basis);
         return Ke;
 }/*}}}*/
@@ -2400 +2408 @@
         IssmDouble  viscosity,newviscosity,oldviscosity;
         IssmDouble  viscosity_overshoot,thickness,Jdet;
-        IssmDouble  D_scalar;
         IssmDouble *xyz_list = NULL;
 
         /*Get problem dimension*/
         element->FindParam(&dim,DomainDimensionEnum);
-        if(dim==2) bsize = 2;
-        else       bsize = 5;
 
         /*Fetch number of nodes and dof for this finite element*/
@@ -2414 +2419 @@
         /*Initialize Element matrix and vectors*/
         ElementMatrix* Ke     = element->NewElementMatrix(HOApproximationEnum);
-        IssmDouble*    B      = xNew<IssmDouble>(bsize*numdof);
-        IssmDouble*    Bprime = xNew<IssmDouble>(bsize*numdof);
-        IssmDouble*    D      = xNewZeroInit<IssmDouble>(bsize*bsize);
+        IssmDouble*    dbasis = xNew<IssmDouble>(dim*numnodes);
 
         /*Retrieve all inputs and parameters*/
@@ -2436 +2439 @@
 
                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-                this->GetBHO(B,element,dim,xyz_list,gauss);
-                this->GetBHOprime(Bprime,element,dim,xyz_list,gauss);
-
+                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
                 element->material->ViscosityHO(&viscosity,dim,xyz_list,gauss,vx_input,vy_input);
                 element->material->ViscosityHO(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input);
 
                 newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
-                D_scalar=2.*newviscosity*gauss->weight*Jdet;
-                for(int i=0;i<bsize;i++) D[i*bsize+i]=D_scalar;
-
-                TripleMultiply(B,bsize,numdof,1,
-                                        D,bsize,bsize,0,
-                                        Bprime,bsize,numdof,0,
-                                        &Ke->values[0],1);
+
+                if(dim==3){
+                        for(int i=0;i<numnodes;i++){
+                                for(int j=0;j<numnodes;j++){
+                                        Ke->values[2*i*2*numnodes+2*j] += gauss->weight*Jdet*newviscosity*(
+                                                                4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[2*numnodes+j]*dbasis[2*numnodes+i]
+                                                                );
+                                        Ke->values[2*i*2*numnodes+2*j+1] += gauss->weight*Jdet*newviscosity*(
+                                                                2.*dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + dbasis[0*numnodes+j]*dbasis[1*numnodes+i]
+                                                                );
+                                        Ke->values[(2*i+1)*2*numnodes+2*j] += gauss->weight*Jdet*newviscosity*(
+                                                                2.*dbasis[0*numnodes+j]*dbasis[1*numnodes+i] + dbasis[1*numnodes+j]*dbasis[0*numnodes+i]
+                                                                );
+                                        Ke->values[(2*i+1)*2*numnodes+2*j+1] += gauss->weight*Jdet*newviscosity*(
+                                                                dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + 4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[2*numnodes+j]*dbasis[2*numnodes+i]
+                                                                );
+                                }
+                        }
+                }
+                else{
+                        for(int i=0;i<numnodes;i++){
+                                for(int j=0;j<numnodes;j++){
+                                        Ke->values[i*numnodes+j] += gauss->weight*Jdet*newviscosity*(
+                                                                4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]
+                                                                );
+                                }
+                        }
+                }
         }
 
@@ -2458 +2480 @@
         delete gauss;
         xDelete<IssmDouble>(xyz_list);
-        xDelete<IssmDouble>(D);
-        xDelete<IssmDouble>(Bprime);
-        xDelete<IssmDouble>(B);
+        xDelete<IssmDouble>(dbasis);
         return Ke;
 }/*}}}*/