Changeset 16889 for issm/trunk-jpl/src/c/analyses/BalancevelocityAnalysis.cpp

Timestamp:

11/22/13 10:36:44 (11 years ago)

Author:

Mathieu Morlighem

Message:

NEW: K matrix bal vel

File:

• issm/trunk-jpl/src/c/analyses/BalancevelocityAnalysis.cpp (modified) (1 diff)

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

@@ -60 +60 @@
 /*Finite Element Analysis*/
 ElementMatrix* BalancevelocityAnalysis::CreateKMatrix(Element* element){/*{{{*/
-        _error_("not implemented yet");
+
+        /*Intermediaries */
+        IssmDouble  dhdt,mb,ms,Jdet;
+        IssmDouble  h,gamma,thickness;
+        IssmDouble  hnx,hny,dhnx[2],dhny[2];
+        IssmDouble *xyz_list = NULL;
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int numnodes = element->GetNumberOfNodes();
+
+        /*Initialize Element matrix and vectors*/
+        ElementMatrix* Ke     = element->NewElementMatrix();
+        IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
+        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
+        IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
+        IssmDouble*    HNx    = xNew<IssmDouble>(numnodes);
+        IssmDouble*    HNy    = xNew<IssmDouble>(numnodes);
+        IssmDouble*    H      = xNew<IssmDouble>(numnodes);
+        IssmDouble*    Nx     = xNew<IssmDouble>(numnodes);
+        IssmDouble*    Ny     = xNew<IssmDouble>(numnodes);
+
+        /*Retrieve all Inputs and parameters: */
+        element->GetVerticesCoordinates(&xyz_list);
+        Input* H_input = element->GetInput(ThicknessEnum); _assert_(H_input);
+        h = element->CharacteristicLength();
+
+        /*Get vector N for all nodes and build HNx and HNy*/
+        element->GetInputListOnNodes(Nx,SurfaceSlopeXEnum);
+        element->GetInputListOnNodes(Ny,SurfaceSlopeYEnum);
+        element->GetInputListOnNodes(H,ThicknessEnum);
+        for(int i=0;i<numnodes;i++){
+                IssmDouble norm=sqrt(Nx[i]*Nx[i]+Ny[i]*Ny[i]+1.e-10);
+                HNx[i] = -H[i]*Nx[i]/norm;
+                HNy[i] = -H[i]*Ny[i]/norm;
+        }
+
+        /*Start looping on the number of gaussian points:*/
+        Gauss* gauss=element->NewGauss(2);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+                gauss->GaussPoint(ig);
+
+                H_input->GetInputValue(&thickness,gauss);
+                if(thickness<50.) thickness=50.;
+                element->ValueP1DerivativesOnGauss(&dhnx[0],HNx,xyz_list,gauss);
+                element->ValueP1DerivativesOnGauss(&dhny[0],HNy,xyz_list,gauss);
+                element->ValueP1OnGauss(&hnx,HNx,gauss);
+                element->ValueP1OnGauss(&hny,HNy,gauss);
+
+                gamma=h/(2.*thickness+1.e-10);
+
+                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                element->NodalFunctions(basis,gauss);
+                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+                for(int i=0;i<numnodes;i++){
+                        for(int j=0;j<numnodes;j++){
+                                Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
+                                                        (basis[i]+gamma*(basis[i]*(dhnx[0]+dhny[1]) + dbasis[0*numnodes+i]*hnx + dbasis[1*numnodes+i]*hny))*
+                                                        (basis[j]*(dhnx[0]+dhny[1])  + dbasis[0*numnodes+j]*hnx + dbasis[1*numnodes+j]*hny)
+                                                        );
+                        }
+                }
+        }
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(xyz_list);
+        xDelete<IssmDouble>(basis);
+        xDelete<IssmDouble>(dbasis);
+        xDelete<IssmDouble>(H);
+        xDelete<IssmDouble>(Nx);
+        xDelete<IssmDouble>(Ny);
+        xDelete<IssmDouble>(HNx);
+        xDelete<IssmDouble>(HNy);
+        xDelete<IssmDouble>(B);
+        delete gauss;
+        return Ke;
 }/*}}}*/
 ElementVector* BalancevelocityAnalysis::CreatePVector(Element* element){/*{{{*/