Changeset 21573

Timestamp:

02/22/17 07:21:50 (8 years ago)

Author:

tsantos

Message:

CHG: changes in FemModel->ReMesh (now it just sets elements, vertices, nodes, constrains and materials). Changes in refpattern file names. NOTE: AMR not working yet.

Location:

issm/trunk-jpl/src/c

Files:

• classes/AdaptiveMeshRefinement.cpp (modified) (1 diff)
• classes/FemModel.cpp (modified) (17 diffs)
• classes/FemModel.h (modified) (2 diffs)
• cores/transient_core.cpp (modified) (1 diff)

issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp

@@ -546 +546 @@
    std::string filename1 = filepath + "/2D_Triang_Rib_3.rpt";
    std::string filename2 = filepath + "/2D_Triang_Rib_4.rpt";
-   std::string filename3 = filepath + "/2D_Triang_Rib2_Side_3_4.rpt";
-   std::string filename4 = filepath + "/2D_Triang_Rib2_Side_3_4permuted.rpt";
-   std::string filename5 = filepath + "/2D_Triang_Rib2_Side_3_5.rpt";
-        std::string filename6 = filepath + "/2D_Triang_Rib2_Side_3_5permuted.rpt";
+   std::string filename3 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_4.rpt";
+   std::string filename4 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_4_permuted.rpt";
+   std::string filename5 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_5.rpt";
+        std::string filename6 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_5_permuted.rpt";
    std::string filename7 = filepath + "/2D_Triang_Rib_5.rpt";
 

issm/trunk-jpl/src/c/classes/FemModel.cpp

@@ -2934 +2934 @@
         /*Get vertices coordinates of the coarse mesh (father mesh)*/
         /*elements comes in Matlab indexing*/
-        this->GetMesh(this,&x,&y,&z,&elements);
+        this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elements);
         
         /*Create initial mesh (coarse mesh) in neopz data structure*/ 
@@ -2953 +2953 @@
 }
 /*}}}*/
-FemModel* FemModel::ReMesh(void){/*{{{*/
+void FemModel::ReMesh(void){/*{{{*/
         
         /*Variables*/
@@ -2973 +2973 @@
         this->ElementsAndVerticesPartitioning(newnumberofvertices,newnumberofelements,elementswidth,newelementslist,&my_elements,&my_vertices);
 
-        /*Creating new femmodel with new mesh*/
-        FemModel* output = new FemModel(*this);
-
-        /*Copy basic attributes:*/
-        output->nummodels        = this->nummodels;
-        output->solution_type    = this->solution_type;
-        output->analysis_counter = this->analysis_counter;
-
-        /*Now, deep copy arrays:*/
-        output->analysis_type_list=xNew<int>(this->nummodels);
-        xMemCpy<int>(output->analysis_type_list,this->analysis_type_list,this->nummodels);
-        output->profiler=static_cast<Profiler*>(this->profiler->copy());
-        output->parameters=static_cast<Parameters*>(this->parameters->Copy());      
         if(this->loads->Size()!=0){
                 _error_("not supported yet");
         }
-        else{
-                output->loads=new Loads();
-        }
-        output->results=new Results();
-        output->amr = NULL;
-        if(IssmComm::GetRank()==0){
-                if(!this->amr) _error_("this->amr is NULL!");
-                output->amr = new AdaptiveMeshRefinement(*this->amr);
-                if(!output->amr) _error_("output->amr is NULL!");
-        }
 
         /*Create vertices*/
-        output->vertices=new Vertices();
-        this->CreateVertices(newnumberofvertices,newnumberofelements,elementswidth,newelementslist,my_vertices,newx,newy,newz,output->vertices);
+        Vertices* new_vertices=new Vertices();
+        this->CreateVertices(newnumberofvertices,newnumberofelements,elementswidth,newelementslist,my_vertices,newx,newy,newz,new_vertices);
  
         /*Creating elements*/
         /*Just Tria in this version*/
-        output->elements=new Elements();
-        this->CreateElements(newnumberofelements,elementswidth,newelementslist,my_elements,output->elements); 
+        Elements* new_elements=new Elements();
+        this->CreateElements(newnumberofelements,elementswidth,newelementslist,my_elements,new_elements);
 
         /*Creating materials*/
-        output->materials=new Materials();
-        this->CreateMaterials(newnumberofelements,my_elements,output->materials);
-        
-        /*Creating nodes*/
+        Materials* new_materials=new Materials();
+        this->CreateMaterials(newnumberofelements,my_elements,new_materials);
+        
+        /*Creating nodes and constraints*/
         /*Just SSA (2D) and P1 in this version*/
-        output->nodes=new Nodes();
-        
+        Nodes* new_nodes=new Nodes();
+        Constraints* new_constraints=new Constraints();
+
         int nodecounter=0;
-        for(int i=0;i<output->nummodels;i++){//create nodes for each analysis in analysis_type_list
-        
-                int analysis_enum = output->analysis_type_list[i];
+        int constraintcounter=0;
+        for(int i=0;i<this->nummodels;i++){//create nodes for each analysis in analysis_type_list
+        
+                int analysis_enum = this->analysis_type_list[i];
                 
                 /*As the domain is 2D, it is not necessary to create nodes for this analysis*/
@@ -3026 +3005 @@
                 if(analysis_enum==StressbalanceVerticalAnalysisEnum) continue;      
                 
-                this->CreateNodes(newnumberofvertices,my_vertices,nodecounter,analysis_enum,output->nodes);
-                this->UpdateElements(newnumberofelements,newelementslist,my_elements,nodecounter,i,output->elements);
-
-                if(output->nodes->Size()) nodecounter = output->nodes->MaximumId();
-        }
-
-        /*Create constraints*/
-        output->constraints=new Constraints();
-        this->CreateConstraints(newnumberofvertices,newnumberofelements,newx,newy,my_vertices,output->constraints);
-        
-        output->elements->Presort();
-        output->nodes->Presort();
-        output->vertices->Presort();
-        output->loads->Presort();
-        output->materials->Presort();
-        output->constraints->Presort();
+                this->CreateNodes(newnumberofvertices,my_vertices,nodecounter,analysis_enum,new_nodes);
+                if(analysis_enum==StressbalanceAnalysisEnum) this->CreateConstraints(newnumberofvertices,newnumberofelements,nodecounter,constraintcounter,newx,newy,my_vertices,new_constraints);
+                this->UpdateElements(newnumberofelements,newelementslist,my_elements,nodecounter,i,new_elements);
+
+
+                if(IssmComm::GetRank()==2){
+                        _printf_("New constraints: CPU 1\n");
+                        //this->constraints->DeepEcho();
+                        new_constraints->DeepEcho();
+                        printarray(my_vertices,newnumberofvertices,1);
+                        _error_("stop!");
+                }
+
+
+                if(new_nodes->Size()) nodecounter=new_nodes->MaximumId();
+                constraintcounter = new_constraints->NumberOfConstraints();
+                /*Make sure nodecounter is at least 0 (if no node exists, maxid will be -1*/
+                _assert_(nodecounter>=0);
+        }
+
+        new_elements->Presort();
+        new_nodes->Presort();
+        new_vertices->Presort();
+        this->loads->Presort();
+        new_materials->Presort();
+        new_constraints->Presort();
 
         /*reset hooks for elements, loads and nodes: */
-        output->elements->ResetHooks();
-        output->loads->ResetHooks();
-        output->materials->ResetHooks();
+        new_elements->ResetHooks();
+        this->loads->ResetHooks();
+        new_materials->ResetHooks();
 
         /*do the post-processing of the datasets to get an FemModel that can actually run analyses: */
         int analysis_type;
-        for(int i=0;i<output->nummodels;i++){
-                analysis_type=output->analysis_type_list[i];
-                output->SetCurrentConfiguration(analysis_type);
-                ConfigureObjectsx(output->elements,output->loads,output->nodes,output->vertices,output->materials,output->parameters);
+        for(int i=0;i<this->nummodels;i++){
+                analysis_type=this->analysis_type_list[i];
+                //SetCurrentConfiguration(analysis_type);
+
+                this->analysis_counter=i;       
+                /*Now, plug analysis_counter and analysis_type inside the parameters: */
+                this->parameters->SetParam(this->analysis_counter,AnalysisCounterEnum);
+                this->parameters->SetParam(analysis_type,AnalysisTypeEnum);
+                this->parameters->SetParam(analysis_type,ConfigurationTypeEnum);
+
+                /*configure elements, loads and nodes, for this new analysis: */
+                new_elements->SetCurrentConfiguration(new_elements,this->loads,new_nodes,new_vertices,new_materials,this->parameters);
+                this->loads->SetCurrentConfiguration(new_elements,this->loads,new_nodes,new_vertices,new_materials,this->parameters);
+
+                /*take care of toolkits options, that depend on this analysis type (present only after model processor)*/
+                if(this->parameters->Exist(ToolkitsOptionsStringsEnum)){
+                        ToolkitsOptionsFromAnalysis(this->parameters,analysis_type);
+                        if(VerboseSolver()) _printf0_("      toolkits Options set for analysis type: " << EnumToStringx(analysis_type) << "\n");
+                }
+                
+                ConfigureObjectsx(new_elements,this->loads,new_nodes,new_vertices,new_materials,this->parameters);
                 if(i==0){ 
-                        VerticesDofx(output->vertices,output->parameters); //only call once, we only have one set of vertices
-                        //GetMaskOfIceVerticesLSMx(output); //itapopo it needs element->inputs
-                }
-                SpcNodesx(output->nodes,output->constraints,output->parameters,analysis_type);
-                NodesDofx(output->nodes,output->parameters,analysis_type);
+                        VerticesDofx(new_vertices,this->parameters); //only call once, we only have one set of vertices
+                }
+                SpcNodesx(new_nodes,new_constraints,this->parameters,analysis_type);
+                NodesDofx(new_nodes,this->parameters,analysis_type);
         }
 
         /*Finally: interpolate all inputs and insert them into the new elements.*/
-        this->InterpolateInputs(output);
-        GetMaskOfIceVerticesLSMx(output);
+        this->InterpolateInputs(new_vertices,new_elements);
+
+        if(IssmComm::GetRank()==2){
+                _printf_("New constraints: CPU 1\n");
+                //this->constraints->DeepEcho();
+                new_constraints->DeepEcho();
+                _error_("stop!");
+        }
+
+
+        /*Delete old structure and set new pointers*/
+        delete this->vertices;          this->vertices          = new_vertices;
+        delete this->elements;          this->elements          = new_elements;
+        delete this->nodes;                     this->nodes                     = new_nodes;
+        delete this->constraints;       this->constraints       = new_constraints;
+        delete this->materials;         this->materials = new_materials;
+        
+        GetMaskOfIceVerticesLSMx(this);
 
         /*Reset current configuration: */
-        analysis_type=output->analysis_type_list[this->analysis_counter];
-        output->SetCurrentConfiguration(analysis_type);
+        analysis_type=this->analysis_type_list[this->analysis_counter];
+        SetCurrentConfiguration(analysis_type);
 
         /*Cleanup*/
@@ -3078 +3099 @@
         xDelete<bool>(my_elements);
 
-        return output;
-
-}
-/*}}}*/
-void FemModel::InterpolateInputs(FemModel* newfemmodel){/*{{{*/
+        return;
+
+}
+/*}}}*/
+void FemModel::InterpolateInputs(Vertices* newfemmodel_vertices,Elements* newfemmodel_elements){/*{{{*/
 
         int maxinputs = MaximumNumberOfDefinitionsEnum;
@@ -3145 +3166 @@
         /*========== Deal with P0 inputs ==========*/
         int numelementsold = this->elements->NumberOfElements();
-        int numelementsnew = newfemmodel->elements->NumberOfElements();
+        int numelementsnew = newfemmodel_elements->NumberOfElements();
         int numverticesold = this->vertices->NumberOfVertices();
-        int numverticesnew = newfemmodel->vertices->NumberOfVertices();
+        int numverticesnew = newfemmodel_vertices->NumberOfVertices();
         IssmDouble* P0inputsold = xNew<IssmDouble>(numelementsold*numP0inputs);
         IssmDouble* P0inputsnew = NULL;
@@ -3186 +3207 @@
         
         /*Get the old mesh*/
-        this->GetMesh(this,&Xold,&Yold,&Zold,&Indexold);
+        this->GetMesh(this->vertices,this->elements,&Xold,&Yold,&Zold,&Indexold);
 
         /*Get the new mesh*/
-        this->GetMesh(newfemmodel,&Xnew,&Ynew,&Znew,&Indexnew);
+        this->GetMesh(newfemmodel_vertices,newfemmodel_elements,&Xnew,&Ynew,&Znew,&Indexnew);
 
         /*Calculate the center points xc and xy*/
@@ -3225 +3246 @@
 
                 /*Update elements from inputs: */
-                for(int j=0;j<newfemmodel->elements->Size();j++){
-                        Element* element=xDynamicCast<Element*>(newfemmodel->elements->GetObjectByOffset(j));
+                for(int j=0;j<newfemmodel_elements->Size();j++){
+                        Element* element=xDynamicCast<Element*>(newfemmodel_elements->GetObjectByOffset(j));
                         switch(P0input_interp[i]){      
                                 case P0Enum:
@@ -3255 +3276 @@
 
                 /*Update elements from inputs: */
-                InputUpdateFromVectorx(newfemmodel,vector,P1input_enums[i],VertexSIdEnum);//VertexSId because vector is serial in SId indexing
-                
+                //InputUpdateFromVectorx(newfemmodel,vector,P1input_enums[i],VertexSIdEnum);//VertexSId because vector is serial in SId indexing
+                for(int j=0;j<newfemmodel_elements->Size();j++){
+                        Element* element=xDynamicCast<Element*>(newfemmodel_elements->GetObjectByOffset(j));
+                        element->InputUpdateFromVector(vector,P1input_enums[i],VertexSIdEnum);
+                }
+
                 xDelete<IssmDouble>(vector);
         }
@@ -3511 +3536 @@
 }
 /*}}}*/
-void FemModel::GetMesh(FemModel* femmodel,IssmDouble** px, IssmDouble** py, IssmDouble** pz, int** pelementslist){/*{{{*/
-
-        if(!femmodel) _error_("GetMesh: NULL femmodel.");
+void FemModel::GetMesh(Vertices* femmodel_vertices, Elements* femmodel_elements,IssmDouble** px, IssmDouble** py, IssmDouble** pz, int** pelementslist){/*{{{*/
+
+        if(!femmodel_vertices) _error_("GetMesh: vertices are NULL.");
+        if(!femmodel_elements) _error_("GetMesh: elements are NULL.");
         
         int numberofvertices, numberofelements;
@@ -3523 +3549 @@
         
         /*Get vertices coordinates*/
-        VertexCoordinatesx(&x, &y, &z, femmodel->vertices,false) ;
-
-        numberofvertices = femmodel->vertices->NumberOfVertices();
-        numberofelements = femmodel->elements->NumberOfElements();
+        VertexCoordinatesx(&x, &y, &z, femmodel_vertices,false) ;
+
+        numberofvertices = femmodel_vertices->NumberOfVertices();
+        numberofelements = femmodel_elements->NumberOfElements();
         
         /*Get element vertices*/
@@ -3535 +3561 @@
 
         /*Go through elements, and for each element, get vertices*/
-   for(int i=0;i<femmodel->elements->Size();i++){
-        Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+   for(int i=0;i<femmodel_elements->Size();i++){
+        Element* element=xDynamicCast<Element*>(femmodel_elements->GetObjectByOffset(i));
         element->GetVerticesSidList(elem_vertices);
         vid1->SetValue(element->sid,elem_vertices[0],INS_VAL);
@@ -3578 +3604 @@
 }
 /*}}}*/
-void FemModel::CreateConstraints(int newnumberofvertices,int newnumberofelements,IssmDouble* newx,IssmDouble* newy,int* my_vertices,Constraints* constraints){/*{{{*/
-
+void FemModel::CreateConstraints(int newnumberofvertices,int newnumberofelements,int nodecounter,int constraintcounter,IssmDouble* newx,IssmDouble* newy,int* my_vertices,Constraints* constraints){/*{{{*/
+
+        /*itapopo ATTENTION: JUST SPCVX AND SPCVY TO TEST!!!*/
+        /*OTHERS CONSTRAINTS MUST BE IMPLEMENTED!!!*/
+        
         /*Get x and y of the mesh i-1*/
         int numberofvertices                    = this->vertices->NumberOfVertices();
@@ -3589 +3618 @@
 
         /*elementslist is in Matlab indexing*/
-        this->GetMesh(this,&x,&y,&z,&elementslist);
-        
-        /*itapopo ATTENTION: JUST SPCVX AND SPCVY TO TEST!!!*/
-        /*OTHERS CONSTRAINTS MUST BE IMPLEMENTED!!!*/
-
+        this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elementslist);
+        
         /*Get spcvx and spcvy for mesh i-1*/
-        IssmDouble *spcvx=NULL;
-        IssmDouble *spcvy=NULL;
-        IssmDouble *spcvxflag=NULL;
-        IssmDouble *spcvyflag=NULL;
-        int numberofnodes_analysistype=this->nodes->NumberOfNodes(StressbalanceAnalysisEnum);
-        Vector<IssmDouble>* vspcvx=new Vector<IssmDouble>(numberofnodes_analysistype);
-        Vector<IssmDouble>* vspcvy=new Vector<IssmDouble>(numberofnodes_analysistype);
-        Vector<IssmDouble>* vspcvxflag=new Vector<IssmDouble>(numberofnodes_analysistype);
-        Vector<IssmDouble>* vspcvyflag=new Vector<IssmDouble>(numberofnodes_analysistype);
+        IssmDouble *spcvx                                               = NULL;
+        IssmDouble *spcvy                                               = NULL;
+        IssmDouble *spcvxflag                           = NULL;
+        IssmDouble *spcvyflag                           = NULL;
+        int numberofnodes_analysistype  = this->nodes->NumberOfNodes(StressbalanceAnalysisEnum);
+        Vector<IssmDouble>* vspcvx                      = new Vector<IssmDouble>(numberofnodes_analysistype);
+        Vector<IssmDouble>* vspcvy                      = new Vector<IssmDouble>(numberofnodes_analysistype);
+        Vector<IssmDouble>* vspcvxflag  = new Vector<IssmDouble>(numberofnodes_analysistype);
+        Vector<IssmDouble>* vspcvyflag  = new Vector<IssmDouble>(numberofnodes_analysistype);
+        
         for(int i=0;i<numberofnodes_analysistype;i++){ 
                 vspcvx->SetValue(i,0.,INS_VAL);
@@ -3612 +3639 @@
 
         for(int i=0;i<this->constraints->Size();i++){
-                SpcStatic* spc=xDynamicCast<SpcStatic*>(this->constraints->GetObjectByOffset(i));
-                int dof=spc->GetDof();
-                int node=spc->GetNodeId();
-                IssmDouble spcvalue=spc->GetValue(); 
-                int nodeindex=node-1;
-                if(dof==0) {//vx
+                SpcStatic* spc                  = xDynamicCast<SpcStatic*>(this->constraints->GetObjectByOffset(i));
+                int dof                                 = spc->GetDof();
+                int node                                        = spc->GetNodeId();
+                IssmDouble spcvalue     = spc->GetValue(); 
+                int nodeindex                   = node-1;
+                
+                if(IssmComm::GetRank()==0)      
+                  {//itapopo
+                 if(dof==0) {//vx
+                        if(IssmComm::GetRank()==0) _printf_("Node: " << nodeindex << "\tdof: " << dof << "\tspcvalue: " << spcvalue << "\n") ;
                         vspcvx->SetValue(nodeindex,spcvalue,INS_VAL);
                         vspcvxflag->SetValue(nodeindex,1,INS_VAL);
                 } 
-                else if(dof==1){//vy
+                if(dof==1){//vy
                         vspcvy->SetValue(nodeindex,spcvalue,INS_VAL);
                         vspcvyflag->SetValue(nodeindex,1,INS_VAL);
                 }
-                else{
-                        /*nothing here*/
-                }
-        }
+
+                  }//itapopo
+        }
+
+        #ifdef _HAVE_PETSC_
+                _printf_("Tem PETSC\n");
+        #endif
 
         /*Assemble*/
@@ -3636 +3670 @@
         vspcvyflag->Assemble();
 
+        IssmDouble pvalue1,pvalue2;
+        vspcvx->GetValue(&pvalue1,0);
+        vspcvxflag->GetValue(&pvalue2,0);
+        if(IssmComm::GetRank()==0) _printf_("Node: " << 0 << "\tspcvalue: " << pvalue1 << "\tflag: " << pvalue2 << "\n") ;
+        vspcvx->GetValue(&pvalue1,5);
+        vspcvxflag->GetValue(&pvalue2,5);
+        if(IssmComm::GetRank()==0) _printf_("Node: " << 5 << "\tspcvalue: " << pvalue1 << "\tflag: " << pvalue2 << "\n") ;
+        
+        
+        
         /*Serialize*/
         spcvx            = vspcvx->ToMPISerial();
@@ -3657 +3701 @@
         InterpFromMeshToMesh2dx(&newspcvxflag,elementslist,x,y,nods_data,nels_data,spcvxflag,M_data,N_data,newx,newy,N_interp,NULL);
         InterpFromMeshToMesh2dx(&newspcvyflag,elementslist,x,y,nods_data,nels_data,spcvyflag,M_data,N_data,newx,newy,N_interp,NULL);
-        
-        int nodecounter                 = 0; //itapopo deve começar pelo primeiro nó do StressbalanceAnalysis
+
+        if(IssmComm::GetRank()==1){
+                _printf_("CPU 1:\n");
+                _printf_("Old spcvx:\n");
+                printarray(spcvx,nods_data,1);
+                _printf_("Old spcvxflag:\n");
+                printarray(spcvxflag,nods_data,1);
+                _error_("Stop!!");
+        }
+
+
+
+
+
         int count                                       = 0;
-        int constraintcounter   = 0; //itapopo
         IssmDouble eps                          = 1.e-8;
 

issm/trunk-jpl/src/c/classes/FemModel.h

@@ -148 +148 @@
                 /*Adaptive mesh refinement methods*/
                 void InitializeAdaptiveRefinement(void);
-                FemModel* ReMesh(void);
-                void GetMesh(FemModel* femmodel,IssmDouble** px, IssmDouble** py, IssmDouble** pz, int** pelementslist);
+                void ReMesh(void);
+                void GetMesh(Vertices* femmodel_vertices,Elements* femmodel_elements,IssmDouble** px, IssmDouble** py, IssmDouble** pz, int** pelementslist);
                 int GetElementsWidth(){return 3;};//just tria elements in this first version
                 void ExecuteRefinement(int &numberofvertices,int &numberofelements,IssmDouble** px,IssmDouble** py,IssmDouble** pz,int** pelementslist);
@@ -157 +157 @@
                 void CreateMaterials(int newnumberofelements,bool* my_elements,Materials* materials);
                 void CreateNodes(int newnumberofvertices,int* my_vertices,int nodecounter,int analysis_enum,Nodes* nodes);
-                void CreateConstraints(int newnumberofvertices,int newnumberofelements,IssmDouble* newx,IssmDouble* newy,int* my_vertices,Constraints* constraints);
-                void InterpolateInputs(FemModel* femmodel);
+                void CreateConstraints(int newnumberofvertices,int newnumberofelements,int nodecounter,int constraintcounter,IssmDouble* newx,IssmDouble* newy,int* my_vertices,Constraints* constraints);
+                void InterpolateInputs(Vertices* newfemmodel_vertices,Elements* newfemmodel_elements);
                 void UpdateElements(int newnumberofelements,int* newelementslist,bool* my_elements,int nodecounter,int analysis_counter,Elements* newelements);
                 void ElementsAndVerticesPartitioning(int& newnumberofvertices,int& newnumberofelements,int& elementswidth,int* newelementslist,bool** pmy_elements,int** pmy_vertices);

issm/trunk-jpl/src/c/cores/transient_core.cpp

@@ -171 +171 @@
                 #ifdef _HAVE_NEOPZ_
                 int isamr=true;
-                if(isamr){
-                        FemModel* newfemmodel = femmodel->ReMesh();
-                        delete femmodel;
-                        femmodel = newfemmodel;
-                }
+                if(isamr) femmodel->ReMesh();
                 #endif