Index: /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp (revision 21572) +++ /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp (revision 21573) @@ -546,8 +546,8 @@ 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"; Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 21572) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 21573) @@ -2934,5 +2934,5 @@ /*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,5 +2953,5 @@ } /*}}}*/ -FemModel* FemModel::ReMesh(void){/*{{{*/ +void FemModel::ReMesh(void){/*{{{*/ /*Variables*/ @@ -2973,52 +2973,31 @@ 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(this->nummodels); - xMemCpy(output->analysis_type_list,this->analysis_type_list,this->nummodels); - output->profiler=static_cast(this->profiler->copy()); - output->parameters=static_cast(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;inummodels;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;inummodels;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,47 +3005,89 @@ 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;inummodels;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;inummodels;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,9 +3099,9 @@ xDelete(my_elements); - return output; - -} -/*}}}*/ -void FemModel::InterpolateInputs(FemModel* newfemmodel){/*{{{*/ + return; + +} +/*}}}*/ +void FemModel::InterpolateInputs(Vertices* newfemmodel_vertices,Elements* newfemmodel_elements){/*{{{*/ int maxinputs = MaximumNumberOfDefinitionsEnum; @@ -3145,7 +3166,7 @@ /*========== 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(numelementsold*numP0inputs); IssmDouble* P0inputsnew = NULL; @@ -3186,8 +3207,8 @@ /*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,6 +3246,6 @@ /*Update elements from inputs: */ - for(int j=0;jelements->Size();j++){ - Element* element=xDynamicCast(newfemmodel->elements->GetObjectByOffset(j)); + for(int j=0;jSize();j++){ + Element* element=xDynamicCast(newfemmodel_elements->GetObjectByOffset(j)); switch(P0input_interp[i]){ case P0Enum: @@ -3255,6 +3276,10 @@ /*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;jSize();j++){ + Element* element=xDynamicCast(newfemmodel_elements->GetObjectByOffset(j)); + element->InputUpdateFromVector(vector,P1input_enums[i],VertexSIdEnum); + } + xDelete(vector); } @@ -3511,7 +3536,8 @@ } /*}}}*/ -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,8 +3549,8 @@ /*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,6 +3561,6 @@ /*Go through elements, and for each element, get vertices*/ - for(int i=0;ielements->Size();i++){ - Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); + for(int i=0;iSize();i++){ + Element* element=xDynamicCast(femmodel_elements->GetObjectByOffset(i)); element->GetVerticesSidList(elem_vertices); vid1->SetValue(element->sid,elem_vertices[0],INS_VAL); @@ -3578,6 +3604,9 @@ } /*}}}*/ -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,19 +3618,17 @@ /*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* vspcvx=new Vector(numberofnodes_analysistype); - Vector* vspcvy=new Vector(numberofnodes_analysistype); - Vector* vspcvxflag=new Vector(numberofnodes_analysistype); - Vector* vspcvyflag=new Vector(numberofnodes_analysistype); + IssmDouble *spcvx = NULL; + IssmDouble *spcvy = NULL; + IssmDouble *spcvxflag = NULL; + IssmDouble *spcvyflag = NULL; + int numberofnodes_analysistype = this->nodes->NumberOfNodes(StressbalanceAnalysisEnum); + Vector* vspcvx = new Vector(numberofnodes_analysistype); + Vector* vspcvy = new Vector(numberofnodes_analysistype); + Vector* vspcvxflag = new Vector(numberofnodes_analysistype); + Vector* vspcvyflag = new Vector(numberofnodes_analysistype); + for(int i=0;iSetValue(i,0.,INS_VAL); @@ -3612,21 +3639,28 @@ for(int i=0;iconstraints->Size();i++){ - SpcStatic* spc=xDynamicCast(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(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,4 +3670,14 @@ 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,8 +3701,19 @@ 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; Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 21572) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 21573) @@ -148,6 +148,6 @@ /*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,6 +157,6 @@ 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); Index: /issm/trunk-jpl/src/c/cores/transient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 21572) +++ /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 21573) @@ -171,9 +171,5 @@ #ifdef _HAVE_NEOPZ_ int isamr=true; - if(isamr){ - FemModel* newfemmodel = femmodel->ReMesh(); - delete femmodel; - femmodel = newfemmodel; - } + if(isamr) femmodel->ReMesh(); #endif