Changeset 16875

Timestamp:

11/21/13 19:51:14 (12 years ago)

Author:

Mathieu Morlighem

Message:

NEW: added KMatrix damage

Location:

issm/trunk-jpl/src/c/analyses

Files:

• DamageEvolutionAnalysis.cpp (modified) (3 diffs)
• DamageEvolutionAnalysis.h (modified) (1 diff)

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

@@ -96 +96 @@
 /*Finite Element Analysis*/
 ElementMatrix* DamageEvolutionAnalysis::CreateKMatrix(Element* element){/*{{{*/
-        _error_("not implemented yet");
-}/*}}}*/
-ElementVector* DamageEvolutionAnalysis::CreatePVector(Element* element){/*{{{*/
 
         /*Intermediaries*/
@@ -118 +115 @@
 
         /*Intermediaries */
+        int         stabilization;
+        IssmDouble  Jdet,dt,D_scalar;
+        IssmDouble  vel,vx,vy,dvxdx,dvydy,dvx[2],dvy[2];
+        IssmDouble *xyz_list  = NULL;
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int numnodes = basalelement->GetNumberOfNodes();
+
+        /*Initialize Element vector*/
+        ElementMatrix* Ke     = basalelement->NewElementMatrix(NoneApproximationEnum);
+        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
+        IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
+        IssmDouble*    Bprime = xNew<IssmDouble>(2*numnodes);
+        IssmDouble     D[2][2]={0.};
+
+        /*Retrieve all inputs and parameters*/
+        basalelement->GetVerticesCoordinates(&xyz_list);
+        basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
+        basalelement->FindParam(&stabilization,DamageStabilizationEnum);
+        Input* vxaverage_input=NULL;
+        Input* vyaverage_input=NULL;
+        if(meshtype==Mesh2DhorizontalEnum){
+                vxaverage_input=basalelement->GetInput(VxEnum); _assert_(vxaverage_input);
+                vyaverage_input=basalelement->GetInput(VyEnum); _assert_(vyaverage_input);
+        }
+        else{
+                vxaverage_input=basalelement->GetInput(VxAverageEnum); _assert_(vxaverage_input);
+                vyaverage_input=basalelement->GetInput(VyAverageEnum); _assert_(vyaverage_input);
+        }
+        IssmDouble h=basalelement->CharacteristicLength();
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=basalelement->NewGauss(2);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+                gauss->GaussPoint(ig);
+
+                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                basalelement->NodalFunctions(basis,gauss);
+                D_scalar=gauss->weight*Jdet;
+
+                vxaverage_input->GetInputValue(&vx,gauss);
+                vyaverage_input->GetInputValue(&vy,gauss);
+                vxaverage_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+                vyaverage_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+
+                TripleMultiply(basis,1,numnodes,1,
+                                        &D_scalar,1,1,0,
+                                        basis,1,numnodes,0,
+                                        &Ke->values[0],1);
+
+                GetB(B,element,xyz_list,gauss);
+                GetBprime(Bprime,element,xyz_list,gauss);
+
+                dvxdx=dvx[0];
+                dvydy=dvy[1];
+                D_scalar=dt*gauss->weight*Jdet;
+
+                D[0][0]=D_scalar*dvxdx;
+                D[1][1]=D_scalar*dvydy;
+                TripleMultiply(B,2,numnodes,1,
+                                        &D[0][0],2,2,0,
+                                        B,2,numnodes,0,
+                                        &Ke->values[0],1);
+
+                D[0][0]=D_scalar*vx;
+                D[1][1]=D_scalar*vy;
+                TripleMultiply(B,2,numnodes,1,
+                                        &D[0][0],2,2,0,
+                                        Bprime,2,numnodes,0,
+                                        &Ke->values[0],1);
+
+                if(stabilization==2){
+                        /*Streamline upwinding*/
+                        vel=sqrt(vx*vx+vy*vy)+1.e-8;
+                        D[0][0]=h/(2*vel)*vx*vx;
+                        D[1][0]=h/(2*vel)*vy*vx;
+                        D[0][1]=h/(2*vel)*vx*vy;
+                        D[1][1]=h/(2*vel)*vy*vy;
+                }
+                else if(stabilization==1){
+                        /*SSA*/
+                        vxaverage_input->GetInputAverage(&vx);
+                        vyaverage_input->GetInputAverage(&vy);
+                        D[0][0]=h/2.0*fabs(vx);
+                        D[0][1]=0.;
+                        D[1][0]=0.;
+                        D[1][1]=h/2.0*fabs(vy);
+                }
+                if(stabilization==1 || stabilization==2){
+                        D[0][0]=D_scalar*D[0][0];
+                        D[1][0]=D_scalar*D[1][0];
+                        D[0][1]=D_scalar*D[0][1];
+                        D[1][1]=D_scalar*D[1][1];
+                        TripleMultiply(Bprime,2,numnodes,1,
+                                                &D[0][0],2,2,0,
+                                                Bprime,2,numnodes,0,
+                                                &Ke->values[0],1);
+                }
+
+        }
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(xyz_list);
+        xDelete<IssmDouble>(basis);
+        xDelete<IssmDouble>(B);
+        xDelete<IssmDouble>(Bprime);
+        delete gauss;
+        if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
+        return Ke;
+}/*}}}*/
+ElementVector* DamageEvolutionAnalysis::CreatePVector(Element* element){/*{{{*/
+
+        /*Intermediaries*/
+        int      meshtype;
+        Element* basalelement;
+
+        /*Get basal element*/
+        element->FindParam(&meshtype,MeshTypeEnum);
+        switch(meshtype){
+                case Mesh2DhorizontalEnum:
+                        basalelement = element;
+                        break;
+                case Mesh3DEnum:
+                        if(!element->IsOnBed()) return NULL;
+                        basalelement = element->SpawnBasalElement();
+                        break;
+                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+        }
+
+        /*Intermediaries */
         IssmDouble  Jdet,dt;
         IssmDouble  f,damage;
@@ -156 +283 @@
         delete gauss;
         return pe;
+}/*}}}*/
+void DamageEvolutionAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+        /*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++){
+                B[numnodes*0+i] = basis[i];
+                B[numnodes*1+i] = basis[i];
+        }
+
+        /*Clean-up*/
+        xDelete<IssmDouble>(basis);
+}/*}}}*/
+void DamageEvolutionAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+        /*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>(2*numnodes);
+        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+        /*Build B': */
+        for(int i=0;i<numnodes;i++){
+                Bprime[numnodes*0+i] = dbasis[0*numnodes+i];
+                Bprime[numnodes*1+i] = dbasis[1*numnodes+i];
+        }
+
+        /*Clean-up*/
+        xDelete<IssmDouble>(dbasis);
+
 }/*}}}*/
 void DamageEvolutionAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/

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

@@ -23 +23 @@
                 ElementMatrix* CreateKMatrix(Element* element);
                 ElementVector* CreatePVector(Element* element);
+                void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
+                void GetBprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
                 void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
                 void InputUpdateFromSolution(IssmDouble* solution,Element* element);