Changeset 22513

Timestamp:

03/09/18 15:23:52 (7 years ago)

Author:

Mathieu Morlighem

Message:

CHG: switching to non matrixmultiply code

File:

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

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

@@ -1307 +1307 @@
         /*Initialize Element matrix and vectors*/
         ElementMatrix* Ke = element->NewElementMatrix(SSAApproximationEnum);
-        IssmDouble*    B  = xNew<IssmDouble>(dim*numdof);
-        IssmDouble*    D  = xNewZeroInit<IssmDouble>(dim*dim);
+        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
 
         /*Retrieve all inputs and parameters*/
@@ -1341 +1340 @@
                 }
 
-                this->GetBSSAFriction(B,element,dim,xyz_list,gauss);
+                element->NodalFunctions(basis,gauss);
                 element->JacobianDeterminant(&Jdet, xyz_list,gauss);
-                for(int i=0;i<dim;i++) D[i*dim+i]=alpha2*gauss->weight*Jdet;
-
-                TripleMultiply(B,dim,numdof,1,
-                                        D,dim,dim,0,
-                                        B,dim,numdof,0,
-                                        &Ke->values[0],1);
+
+                if(dim==2){
+                        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];
+                                }
+                        }
+                }
         }
 
@@ -1358 +1367 @@
         delete friction;
         xDelete<IssmDouble>(xyz_list);
-        xDelete<IssmDouble>(B);
-        xDelete<IssmDouble>(D);
+        xDelete<IssmDouble>(basis);
         return Ke;
 }/*}}}*/
@@ -1434 +1442 @@
 
         /*Intermediaries*/
-        int         dim,domaintype,bsize;
+        int         dim,domaintype;
         IssmDouble  viscosity,newviscosity,oldviscosity;
         IssmDouble  viscosity_overshoot,thickness,Jdet;
-        IssmDouble  D_scalar;
         IssmDouble *xyz_list = NULL;
 
@@ -1443 +1450 @@
         element->FindParam(&domaintype,DomainTypeEnum);
         switch(domaintype){
-                case Domain2DverticalEnum:   dim = 1; bsize = 1; break;
-                case Domain2DhorizontalEnum: dim = 2; bsize = 3; break;
-                case Domain3DEnum:           dim = 2; bsize = 3; break;
+                case Domain2DverticalEnum:   dim = 1; break;
+                case Domain2DhorizontalEnum: dim = 2; break;
+                case Domain3DEnum:           dim = 2; break;
                 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
         }
@@ -1455 +1462 @@
         /*Initialize Element matrix and vectors*/
         ElementMatrix* Ke     = element->NewElementMatrix(SSAApproximationEnum);
-        IssmDouble*    B      = xNew<IssmDouble>(bsize*numdof);
-        IssmDouble*    Bprime = xNew<IssmDouble>(bsize*numdof);
-        IssmDouble*    D      = xNewZeroInit<IssmDouble>(bsize*bsize);
+        IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
 
         /*Retrieve all inputs and parameters*/
@@ -1478 +1483 @@
 
                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-                this->GetBSSA(B,element,dim,xyz_list,gauss);
-                this->GetBSSAprime(Bprime,element,dim,xyz_list,gauss);
-
+                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+                thickness_input->GetInputValue(&thickness, gauss);
                 element->material->ViscositySSA(&viscosity,dim,xyz_list,gauss,vx_input,vy_input);
                 element->material->ViscositySSA(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input);
-                thickness_input->GetInputValue(&thickness, gauss);
-
                 newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
-                D_scalar=2.*newviscosity*thickness*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==2){
+                        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*thickness*(
+                                                                4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]
+                                                                );
+                                        Ke->values[2*i*2*numnodes+2*j+1] += gauss->weight*Jdet*newviscosity*thickness*(
+                                                                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*thickness*(
+                                                                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*thickness*(
+                                                                dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + 4.*dbasis[1*numnodes+j]*dbasis[1*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*thickness*(
+                                                                4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i]
+                                                                );
+                                }
+                        }
+                }
         }
 
@@ -1501 +1525 @@
         delete gauss;
         xDelete<IssmDouble>(xyz_list);
-        xDelete<IssmDouble>(D);
-        xDelete<IssmDouble>(Bprime);
-        xDelete<IssmDouble>(B);
+        xDelete<IssmDouble>(dbasis);
         return Ke;
 }/*}}}*/