Changeset 21524

Timestamp:

02/06/17 10:39:11 (8 years ago)

Author:

tsantos

Message:

CHG: changes in AMR methods: added new method to get mesh by femmodel; changes in InterpolateInputs (NOTE: AMR not working yet).

Location:

issm/trunk-jpl/src/c/classes

Files:

• FemModel.cpp (modified) (5 diffs)
• FemModel.h (modified) (1 diff)

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

@@ -3103 +3103 @@
         output->materials->ResetHooks();
 
-        printf("-------------- file: FemModel.cpp line: %i\n",__LINE__); 
-        /*Finally: interpolate all inputs*/
-        this->InterpolateInputs(output);
-
-        printf("-------------- file: FemModel.cpp line: %i\n",__LINE__); 
         /*do the post-processing of the datasets to get an FemModel that can actually run analyses: */
         int analysis_type;
@@ -3122 +3117 @@
         }
 
+        /*Finally: interpolate all inputs*/
+        printf("-------------- file: FemModel.cpp line: %i\n",__LINE__); 
+        this->InterpolateInputs(output);
+        
         /*Reset current configuration: */
         analysis_type=output->analysis_type_list[this->analysis_counter];
@@ -3214 +3213 @@
         printarray(P0inputsold,numelementsold,numP0inputs);
 
-        _error_("stop");
+        /*========== Deal with P1 inputs ==========*/
+        IssmDouble* P1inputsold = xNew<IssmDouble>(numverticesold*numP1inputs);
+        IssmDouble* P1inputsnew = NULL;
+        vector      = NULL;
+
+        for(int i=0;i<numP1inputs;i++){
+                printf("Dealing with %s (type: %s)\n",EnumToStringx(P1input_enums[i]),EnumToStringx(P1input_interp[i]));
+                GetVectorFromInputsx(&vector,this,P1input_enums[i],VertexSIdEnum);
+
+                /*Copy vector in matrix*/
+                for(int j=0;j<numverticesold;j++) P1inputsold[j*numP1inputs+i] = vector[j];
+                xDelete<IssmDouble>(vector);
+        }
+        printf("Printing array to make sure it looks alright\n");
+        printarray(P1inputsold,numverticesold,numP1inputs);
+        
         /*Old mesh coordinates*/
         IssmDouble *Xold     = NULL;
         IssmDouble *Yold     = NULL;
+        IssmDouble *Zold                = NULL;
         int        *Indexold = NULL;
+        int        *Indexnew = NULL;
         IssmDouble *Xnew     = NULL;
         IssmDouble *Ynew     = NULL;
-
-
+        IssmDouble *Znew                = NULL;
+        
+        /*Get the old mesh*/
         printf("-------------- file: FemModel.cpp line: %i\n",__LINE__); 
+        this->GetMesh(this,&Xold,&Yold,&Zold,&Indexold);
+
+        /*Get the new mesh*/
+        printf("-------------- file: FemModel.cpp line: %i\n",__LINE__); 
+        this->GetMesh(newfemmodel,&Xnew,&Ynew,&Znew,&Indexnew);
+        
+        printf("-------------- file: FemModel.cpp line: %i\n",__LINE__); 
+        /*Interplate P0 inputs in the new mesh*/
         InterpFromMeshToMesh2dx(&P0inputsnew,Indexold,Xold,Yold,numverticesold,numelementsold,
                                 P0inputsold,numelementsold,numP0inputs,
                                 Xnew,Ynew,numelementsnew,NULL);
-
+        
+        printf("Printing array AFTER INTERP - P0 INPUTS - to make sure it looks alright\n");
+        printarray(P0inputsnew,numelementsnew,numP0inputs);
+        
+        /*Interpolate P1 inputs in the new mesh*/
+        InterpFromMeshToMesh2dx(&P1inputsnew,Indexold,Xold,Yold,numverticesold,numelementsold,
+                                P1inputsold,numverticesold,numP1inputs,
+                                Xnew,Ynew,numverticesnew,NULL);
+
+        printf("Printing array AFTER INTERP - P1 INPUTS - to make sure it looks alright\n");
+        printarray(P1inputsnew,numverticesnew,numP1inputs);
+        
         _error_("STOP");
 
         xDelete<IssmDouble>(P0inputsold);
         xDelete<IssmDouble>(P0inputsnew);
-
+        xDelete<IssmDouble>(P1inputsold);
+        xDelete<IssmDouble>(P1inputsnew);
 
         _error_("STOP");
@@ -3411 +3448 @@
 
         int lid=0;
-
+//itapopo use the GetMesh!!
         for(int j=0;j<newnumberofvertices;j++){
                 if(my_vertices[j]){                             
@@ -3455 +3492 @@
 }
 /*}}}*/
+void FemModel::GetMesh(FemModel* femmodel,IssmDouble** px, IssmDouble** py, IssmDouble** pz, int** pelementslist){/*{{{*/
+
+        if(!femmodel) _error_("GetMesh: NULL femmodel.");
+        
+        int numberofvertices, numberofelements;
+        int elementswidth = 3; //itapopo just 2D mesh in this version (just tria elements)
+        IssmDouble *x           = NULL;
+        IssmDouble *y           = NULL;
+        IssmDouble *z           = NULL; 
+        int* elementslist = NULL;
+        
+        /*Get vertices coordinates*/
+        VertexCoordinatesx(&x, &y, &z, femmodel->vertices,false) ;
+
+        numberofvertices = femmodel->vertices->NumberOfVertices();
+        numberofelements = femmodel->elements->NumberOfElements();
+        
+        /*Get element vertices*/
+        int* elem_vertices=xNew<int>(elementswidth);
+        Vector<IssmDouble>* vid1= new Vector<IssmDouble>(numberofelements);
+        Vector<IssmDouble>* vid2= new Vector<IssmDouble>(numberofelements);
+        Vector<IssmDouble>* vid3= new Vector<IssmDouble>(numberofelements);
+
+        /*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));
+        element->GetVerticesSidList(elem_vertices);
+        vid1->SetValue(element->sid,elem_vertices[0],INS_VAL);
+        vid2->SetValue(element->sid,elem_vertices[1],INS_VAL);
+        vid3->SetValue(element->sid,elem_vertices[2],INS_VAL);
+   }
+                
+        /*Assemble*/
+   vid1->Assemble();
+   vid2->Assemble();
+   vid3->Assemble();
+
+   /*Serialize*/
+        IssmDouble *id1 = vid1->ToMPISerial();
+   IssmDouble *id2 = vid2->ToMPISerial();
+        IssmDouble *id3 = vid3->ToMPISerial();
+        
+        /*Construct elements list*/
+        elementslist=xNew<int>(numberofelements*elementswidth);
+        if(numberofelements*elementswidth<0) _error_("numberofelements negative.");
+        for(int i=0;i<numberofelements;i++){
+                std::cout << (int)id1[i]+1 << " " << (int)id2[i]+1 <<" " << (int)id3[i]+1 <<  std::endl;
+                elementslist[elementswidth*i+0] = (int)id1[i]+1; //InterpMesh wants Matlab indexing
+                elementslist[elementswidth*i+1] = (int)id2[i]+1; //InterpMesh wants Matlab indexing
+                elementslist[elementswidth*i+2] = (int)id3[i]+1; //InterpMesh wants Matlab indexinf
+        }
+        
+        /*Assign pointers*/
+        *px                             = x;
+        *py                             = y;
+        *pz                             = z;
+        *pelementslist = elementslist;
+}
+/*}}}*/
 void FemModel::CreateConstraints(int newnumberofvertices,int newnumberofelements,int* newelementslist,IssmDouble* newx,IssmDouble* newy,int* my_vertices,Constraints* constraints){/*{{{*/
 

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

@@ -149 +149 @@
                 void InitializeAdaptiveRefinement(void);
                 FemModel* ReMesh(void);
+                void GetMesh(FemModel* femmodel,IssmDouble** px, IssmDouble** py, IssmDouble** pz, int** pelementslist);
                 void ExecuteRefinement(int &newnumberofvertices,int &newnumberofelements,IssmDouble** newx,IssmDouble** newy,IssmDouble** newz,int** newelementslist);
                 void GetMaskLevelSet(IssmDouble** pmasklevelset);