Changeset 17334

Timestamp:

02/21/14 10:24:59 (11 years ago)

Author:

jbondzio

Message:

ADD: LsfReinitializationAnalysis: CreateKMatrix, CreatePVector. Preparing Picard Iteration on Levelsetfunction for initialization.

Location:

issm/trunk-jpl/src

Files:

• c/analyses/LsfReinitializationAnalysis.cpp (modified) (7 diffs)
• c/analyses/LsfReinitializationAnalysis.h (modified) (1 diff)
• c/classes/Elements/Element.cpp (modified) (1 diff)
• c/shared/Enum/EnumDefinitions.h (modified) (1 diff)
• c/shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• c/shared/Enum/StringToEnumx.cpp (modified) (1 diff)
• m/enum/EnumDefinitions.py (modified) (1 diff)
• m/enum/LevelsetfunctionPicardEnum.m (added)

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

@@ -5 +5 @@
 #include "../modules/modules.h"
 
+#include "../modules/GetVectorFromInputsx/GetVectorFromInputsx.h"
+
 /*Model processing*/
 int  LsfReinitializationAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
@@ -13 +15 @@
 }/*}}}*/
 void LsfReinitializationAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
-        /* Do nothing for now */
+        int    finiteelement;
+
+        /*Finite element type*/
+        finiteelement = P1Enum;
+
+        /*Update elements: */
+        int counter=0;
+        for(int i=0;i<iomodel->numberofelements;i++){
+                if(iomodel->my_elements[i]){
+                        Element* element=(Element*)elements->GetObjectByOffset(counter);
+                        element->Update(i,iomodel,analysis_counter,analysis_type,finiteelement);
+                        counter++;
+                }
+        }
 }/*}}}*/
 void LsfReinitializationAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
-        /* Do nothing for now */
+        int finiteelement=P1Enum;
+        ::CreateNodes(nodes,iomodel,LsfReinitializationAnalysisEnum,finiteelement);
 }/*}}}*/
 void LsfReinitializationAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
@@ -26 +42 @@
 
 /*Finite element Analysis*/
-void           LsfReinitializationAnalysis::Core(FemModel* femmodel){/*{{{*/
-        _error_("not implemented yet");
+void  LsfReinitializationAnalysis::Core(FemModel* femmodel){/*{{{*/
+
+        /*parameters: */
+        bool save_results;
+        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
+
+        /*activate formulation: */
+        femmodel->SetCurrentConfiguration(LsfReinitializationAnalysisEnum);
+
+        /* set spcs for reinitialization */
+        if(VerboseSolution()) _printf0_("Update spcs for reinitialization:\n");
+        UpdateReinitSPCs(femmodel);
+
+        if(VerboseSolution()) _printf0_("call computational core for reinitialization:\n");
+//      solutionsequence_lsfreinit_linear(femmodel);
+
+        if(save_results){
+                if(VerboseSolution()) _printf0_("   saving results\n");
+                int outputs = MaskIceLevelsetEnum;
+                femmodel->RequestedOutputsx(&femmodel->results,&outputs,1);
+        }
+
 }/*}}}*/
 ElementVector* LsfReinitializationAnalysis::CreateDVector(Element* element){/*{{{*/
-        _error_("not implemented yet");
+        /*Default, return NULL*/
+        return NULL;
 }/*}}}*/
 ElementMatrix* LsfReinitializationAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
@@ -36 +73 @@
 }/*}}}*/
 ElementMatrix* LsfReinitializationAnalysis::CreateKMatrix(Element* element){/*{{{*/
-        _error_("not implemented yet");
+
+        /*Intermediaries */
+        const int dim = 2;
+        int        i,row,col,stabilization;
+        IssmDouble Jdet,D_scalar,h;
+        IssmDouble dlsf[3],normal[3];
+        IssmDouble norm_dlsf;
+        IssmDouble hx,hy,hz,kappa;
+        IssmDouble* xyz_list = NULL;
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int numnodes = element->GetNumberOfNodes();
+
+        /*Initialize Element vector and other vectors*/
+        ElementMatrix* Ke     = element->NewElementMatrix();
+        IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
+        IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
+        IssmDouble     D[dim][dim];
+
+        /*Retrieve all inputs and parameters*/
+        Input* lsfpicard_input=element->GetInput(LevelsetfunctionPicardEnum); _assert_(lsfpicard_input);
+        element->GetVerticesCoordinates(&xyz_list);
+        h = element->CharacteristicLength();
+
+        /* 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);
+
+                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                GetB(B,element,xyz_list,gauss);
+                GetBprime(Bprime,element,xyz_list,gauss);
+
+                /* Get normal from last iteration on lsf */
+                lsfpicard_input->GetInputDerivativeValue(&dlsf[0],xyz_list,gauss);
+
+                norm_dlsf=0.;
+                for(i=0;i<dim;i++) norm_dlsf+=dlsf[i]*dlsf[i]; 
+                norm_dlsf=sqrt(norm_dlsf); _assert_(norm_dlsf>0.);
+                for(i=0;i<dim;i++)
+                        normal[i]=dlsf[i]/norm_dlsf;
+                
+                D_scalar=gauss->weight*Jdet;
+
+                for(row=0;row<dim;row++)
+                        for(col=0;col<dim;col++)
+                                if(row==col)
+                                        D[row][col]=D_scalar*normal[row];
+                                else
+                                        D[row][col]=0.;
+                TripleMultiply(B,dim,numnodes,1,
+                                        &D[0][0],dim,dim,0,
+                                        Bprime,dim,numnodes,0,
+                                        &Ke->values[0],1);
+
+                /* Stabilization *//*{{{*/
+                stabilization=1;
+                if (stabilization==0){/* no stabilization, do nothing*/}
+                else if(stabilization==1){
+                        /* Artificial Diffusion */
+                        element->ElementSizes(&hx,&hy,&hz);
+                        h=sqrt( pow(hx*normal[0],2) + pow(hy*normal[1],2));
+                        kappa=h/2.; 
+                        D[0][0]=D_scalar*kappa;
+                        D[0][1]=0.;
+                        D[1][0]=0.;
+                        D[1][1]=D_scalar*kappa;
+                        TripleMultiply(Bprime,dim,numnodes,1,
+                                                &D[0][0],dim,dim,0,
+                                                Bprime,dim,numnodes,0,
+                                                &Ke->values[0],1);
+                }
+                else if(stabilization==2){
+                        /*Streamline upwinding - do not use this for extrapolation: yields oscillating results due to smoothing along normal, not across */
+                        for(row=0;row<dim;row++)
+                                for(col=0;col<dim;col++)
+                                        D[row][col]=h/(2.*1.)*normal[row]*normal[col];
+
+                        TripleMultiply(Bprime,dim,numnodes,1,
+                                                &D[0][0],dim,dim,0,
+                                                Bprime,dim,numnodes,0,
+                                                &Ke->values[0],1);
+                }/*}}}*/
+        }/*}}}*/
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(xyz_list);
+        xDelete<IssmDouble>(B);
+        xDelete<IssmDouble>(Bprime);
+        delete gauss;
+        return Ke;
 }/*}}}*/
 ElementVector* LsfReinitializationAnalysis::CreatePVector(Element* element){/*{{{*/
-        _error_("not implemented yet");
-}/*}}}*/
+        
+        /*Intermediaries */
+        IssmDouble Jdet;
+        IssmDouble* xyz_list = NULL;
+        
+        /*Fetch number of nodes */
+        int numnodes = element->GetNumberOfNodes();
+
+        IssmDouble* basis = xNew<IssmDouble>(numnodes);
+        element->GetVerticesCoordinates(&xyz_list);
+
+        /*Initialize Element vector*/
+        ElementVector* pe = element->NewElementVector();
+
+        Gauss* gauss=element->NewGauss(2);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+                gauss->GaussPoint(ig);
+
+                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                element->NodalFunctions(basis,gauss);
+
+                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*basis[i]; 
+        }
+
+        xDelete<IssmDouble>(basis);
+        xDelete<IssmDouble>(xyz_list);
+        return pe;
+        }/*}}}*/
 void LsfReinitializationAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
         _error_("not implemented yet");
@@ -50 +203 @@
         _error_("not implemented yet");
 }/*}}}*/
+void LsfReinitializationAnalysis::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 LsfReinitializationAnalysis::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);
+
+}/*}}}*/
 
 /* Other */
+void LsfReinitializationAnalysis::UpdateReinitSPCs(FemModel* femmodel){/*{{{*/
+
+        int i,k, numnodes;
+        Element* element;
+        Node* node;
+
+        /* deactivate all spcs */
+        for(i=0;i<femmodel->elements->Size();i++){
+                element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+                for(k=0;k<element->GetNumberOfNodes();k++){
+                                node=element->GetNode(k);
+                                node->DofInFSet(0); 
+                }
+        }
+
+        SetDistanceOnIntersectedElements(femmodel);
+
+        /* reactivate spcs on elements intersected by zero levelset */
+        for(i=0;i<femmodel->elements->Size();i++){
+                element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+                if(element->IsZeroLevelset(MaskIceLevelsetEnum)){
+                        /*iterate over nodes and set spc */
+                        numnodes=element->GetNumberOfNodes();
+                        IssmDouble* lsf = xNew<IssmDouble>(numnodes);
+                        element->GetInputListOnNodes(&lsf[0],MaskIceLevelsetEnum);
+                        for(k=0;k<numnodes;k++){
+                                node=element->GetNode(k);
+                                node->ApplyConstraint(1,lsf[k]);
+                        }
+                        xDelete<IssmDouble>(lsf);
+                }
+        }
+
+}/*}}}*/
 void LsfReinitializationAnalysis::SetDistanceOnIntersectedElements(FemModel* femmodel){/*{{{*/
 
         /* Intermediaries */
-        int i;
+        int i,k;
 
         /*Initialize vector with number of vertices*/
         int numvertices=femmodel->vertices->NumberOfVertices();
-        Vector<IssmDouble>* vec_dist_zerolevelset=new Vector<IssmDouble>(numvertices); //vertices that have ice at next time step
-
-        /*TODO: Fill vector with values of old level set function: */
-        for(i=0;i<numvertices;i++){
-                vec_dist_zerolevelset->SetValue(i,0.,INS_VAL);
-        }
+        Element* element;
+
+        Vector<IssmDouble>* vec_dist_zerolevelset = NULL;
+        GetVectorFromInputsx(&vec_dist_zerolevelset, femmodel, MaskIceLevelsetEnum, VertexEnum);
+        
         for(i=0;i<femmodel->elements->Size();i++){
-                Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+                element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+                if(element->IsZeroLevelset(MaskIceLevelsetEnum))
+                        for(k=0;k<element->GetNumberOfVertices();k++)
+                                vec_dist_zerolevelset->SetValue(element->vertices[k]->Sid(),NAN,INS_VAL); 
+        }
+
+        for(i=0;i<femmodel->elements->Size();i++){
+                element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
                 if(element->IsZeroLevelset(MaskIceLevelsetEnum))
                         SetDistanceToZeroLevelsetElement(vec_dist_zerolevelset, element);
@@ -103 +351 @@
 
         /* get sign of levelset function */
-        for(i=0;i<numvertices;i++){
-                if(lsf[i]==0.)
-                        sign_lsf[i]=1.;
-                else
-                        sign_lsf[i]=lsf[i]/fabs(lsf[i]);
-        }
+        for(i=0;i<numvertices;i++)
+                sign_lsf[i]=(lsf[i]>=0.?1.:-1.);
 
         element->ZeroLevelsetCoordinates(&xyz_list_zero, xyz_list, MaskIceLevelsetEnum);
@@ -127 +371 @@
         for(i=0;i<numvertices;i++){
                 vec_signed_dist->GetValue(&lsf_old, element->vertices[i]->Sid());
-                if(fabs(signed_dist[i])<fabs(lsf_old))
+                if(isnan(lsf_old) || fabs(signed_dist[i])<fabs(lsf_old))
                         vec_signed_dist->SetValue(element->vertices[i]->Sid(),signed_dist[i],INS_VAL);
         }

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

@@ -29 +29 @@
         void InputUpdateFromSolution(IssmDouble* solution,Element* element);
         void UpdateConstraints(FemModel* femmodel);
+        void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
+        void GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
         
         /* Other */
+        void UpdateReinitSPCs(FemModel* femmodel);
         void SetDistanceOnIntersectedElements(FemModel* femmodel);
         void SetDistanceToZeroLevelsetElement(Vector<IssmDouble>* vec_dist, Element* element);

issm/trunk-jpl/src/c/classes/Elements/Element.cpp

@@ -434 +434 @@
                                 name==LevelsetfunctionSlopeXEnum ||
                                 name==LevelsetfunctionSlopeYEnum ||
+                                name==LevelsetfunctionPicardEnum ||
                                 name==GradientEnum ||
                                 name==OldGradientEnum  ||

issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

@@ -682 +682 @@
         LevelsetfunctionSlopeXEnum,
         LevelsetfunctionSlopeYEnum,
+        LevelsetfunctionPicardEnum,
         /*}}}*/
         MaximumNumberOfDefinitionsEnum

issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

@@ -641 +641 @@
                 case LevelsetfunctionSlopeXEnum : return "LevelsetfunctionSlopeX";
                 case LevelsetfunctionSlopeYEnum : return "LevelsetfunctionSlopeY";
+                case LevelsetfunctionPicardEnum : return "LevelsetfunctionPicard";
                 case MaximumNumberOfDefinitionsEnum : return "MaximumNumberOfDefinitions";
                 default : return "unknown";

issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

@@ -656 +656 @@
               else if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum;
               else if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum;
+              else if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum;
               else if (strcmp(name,"MaximumNumberOfDefinitions")==0) return MaximumNumberOfDefinitionsEnum;
          else stage=7;

issm/trunk-jpl/src/m/enum/EnumDefinitions.py

@@ -633 +633 @@
 def LevelsetfunctionSlopeXEnum(): return StringToEnum("LevelsetfunctionSlopeX")[0]
 def LevelsetfunctionSlopeYEnum(): return StringToEnum("LevelsetfunctionSlopeY")[0]
+def LevelsetfunctionPicardEnum(): return StringToEnum("LevelsetfunctionPicard")[0]
 def MaximumNumberOfDefinitionsEnum(): return StringToEnum("MaximumNumberOfDefinitions")[0]