Changeset 25611

Timestamp:

09/29/20 11:26:14 (5 years ago)

Author:

Mathieu Morlighem

Message:

CHG: working on MLHO

File:

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

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

@@ -2728 +2728 @@
         /*Fetch number of nodes and dof for this finite element*/
         int numnodes = basalelement->GetNumberOfNodes();
-        int numdof   = numnodes*2;
 
         /*Initialize Element matrix and vectors*/
         ElementMatrix* Ke     = basalelement->NewElementMatrix(MLHOApproximationEnum);
-        IssmDouble*    B      = xNew<IssmDouble>(3*numdof);
-        IssmDouble*    Bprime = xNew<IssmDouble>(3*numdof);
-        IssmDouble*    D      = xNewZeroInit<IssmDouble>(3*3);
+        IssmDouble*    dbasis = xNew<IssmDouble>(3*numnodes);
 
         /*Retrieve all inputs and parameters*/
@@ -2749 +2746 @@
 
                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-                this->GetBSSA(B,basalelement,2,xyz_list,gauss_base);
-                this->GetBSSAprime(Bprime,basalelement,2,xyz_list,gauss_base);
+                basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss_base);
 
                 element->material->ViscosityL1L2(&viscosity,xyz_list,gauss,vx_input,vy_input,surface_input);
 
-                for(int i=0;i<3;i++) D[i*3+i]=2*viscosity*gauss->weight*Jdet;
-
-                TripleMultiply(B,3,numdof,1,
-                                        D,3,3,0,
-                                        Bprime,3,numdof,0,
-                                        &Ke->values[0],1);
+                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*viscosity*(
+                                                        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*viscosity*(
+                                                        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*viscosity*(
+                                                        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*viscosity*(
+                                                        dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + 4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i]
+                                                        );
+                        }
+                }
         }
 
@@ -2770 +2776 @@
         basalelement->DeleteMaterials(); delete basalelement;
         xDelete<IssmDouble>(xyz_list);
-        xDelete<IssmDouble>(D);
-        xDelete<IssmDouble>(Bprime);
-        xDelete<IssmDouble>(B);
+        xDelete<IssmDouble>(dbasis);
         return Ke;
 }/*}}}*/