Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 21515) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 21516) @@ -264,5 +264,6 @@ ./modules/ConstraintsStatex/RiftConstraintsState.cpp\ ./modules/ModelProcessorx/CreateOutputDefinitions.cpp\ - ./modules/OutputDefinitionsResponsex/OutputDefinitionsResponsex.cpp\ + ./modules/OutputDefinitionsResponsex/OutputDefinitionsResponsex.cpp\ + ./modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp\ ./classes/Inputs/PentaInput.cpp\ ./classes/Inputs/TetraInput.cpp @@ -549,5 +550,4 @@ ./modules/InterpFromMesh2dx/InterpFromMesh2dx.cpp\ ./modules/InterpFromMesh2dx/InterpFromMesh2dxt.cpp\ - ./modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp\ ./modules/InterpFromMeshToMesh3dx/InterpFromMeshToMesh3dx.cpp\ ./modules/InterpFromMeshToGridx/InterpFromMeshToGridx.cpp\ Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 21515) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 21516) @@ -84,9 +84,7 @@ this->InitializeAdaptiveRefinement(); FemModel *Test=this->ReMesh();//itapopo: just to test!; - printf(" AFTER REMESH!!!\n"); + //Test->elements->DeepEcho(); //Test->CleanUp(); - printf(" AFTER CLEANUP!!!\n"); //delete Test; - printf(" AFTER DELETE!!!\n"); #endif @@ -3012,81 +3010,50 @@ FemModel* FemModel::ReMesh(void){/*{{{*/ - /*All indexing here is in C type: 0..n-1*/ - int my_rank=IssmComm::GetRank(); - const int elementswidth=3;//just 2D mesh, tria elements - - /*Solutions which will be used to refine the elements*/ - IssmDouble *vx=NULL; //This will be used in constraints - IssmDouble *vy=NULL; //This will be used in constraints - IssmDouble *masklevelset=NULL; - - this->GetMaskLevelSet(&masklevelset); - // #1: TEST MASK LEVEL SET - _printf_(" PRINTING MASKLEVELSET... \n"); - int numberofvertices=this->vertices->NumberOfVertices(); - if(my_rank==0) for(int i=0;iamr->ExecuteRefinement(type_process,vx,vy,masklevelset,newnumberofvertices,newnumberofelements,newnumberofsegments,&newx,&newy,&newz,&newelements,&newsegments); - - if(newnumberofvertices<=0 || newnumberofelements<=0 /*|| newnumberofsegments<=0*/) _error_("Error in the refinement process."); - - /*Cleanup masklevetset*/ - xDelete(masklevelset); - - // #2: TEST NEW MESH - if(my_rank==0){ - _printf_(" PRINTING COORDINATES... \n"); - for(int i=0;i(newnumberofelements*elementswidth); - for(int i=0;i1) _error_("Multithreading not tested yet. Run with 1 cpu."); + + /*Variables*/ + IssmDouble *newx=NULL; + IssmDouble *newy=NULL; + IssmDouble *newz=NULL; + int *newelementslist=NULL; + int newnumberofvertices=-1; + int newnumberofelements=-1; bool* my_elements=NULL; int* my_vertices=NULL; - + int elementswidth=3;//just tria elements in this version + + /*Execute refinement and get the new mesh*/ + this->ExecuteRefinement(newnumberofvertices,newnumberofelements,&newx,&newy,&newz,&newelementslist); + + /*Partitioning the new mesh. Maybe ElementsAndVerticesPartitioning.cpp could be modified to set this without iomodel.*/ this->ElementsAndVerticesPartitioning(newnumberofvertices,newnumberofelements,elementswidth,newelementslist,&my_elements,&my_vertices); /*Creating new femmodel with new mesh*/ - FemModel* output=NULL; - output=new FemModel(*this); + FemModel* output=new FemModel(*this); /*Copy basic attributes:*/ - output->nummodels=this->nummodels; - output->solution_type=this->solution_type; - output->analysis_counter=this->analysis_counter; + 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(); + } + + //itapopo + // output->results=static_cast(this->results->Copy()); + /*Create vertices*/ output->vertices=new Vertices(); this->CreateVertices(newnumberofvertices,newnumberofelements,elementswidth,newelementslist,my_vertices,newx,newy,newz,output->vertices); - - /* #4: TEST THE VERTICES!*/ - //_printf_(" Old vertices deep echo: \n"); - //this->vertices->DeepEcho(); - //_printf_(" New vertices deep echo: \n"); - //output->vertices->DeepEcho(); /*Creating elements*/ @@ -3094,211 +3061,211 @@ output->elements=new Elements(); this->CreateElements(newnumberofelements,elementswidth,newelementslist,my_elements,output->elements); - /*Cleanup*/ - xDelete(newelementslist); +// xDelete(newelementslist); itapopo /*Creating materials*/ output->materials=new Materials(); this->CreateMaterials(newnumberofelements,my_elements,output->materials); + /*Creating nodes*/ /*Just SSA (2D) and P1 in this version*/ output->nodes=new Nodes(); - + int nodecounter=0; - int lid=0; + + for(int i=0;inummodels;i++){//create nodes for each analysis in analysis_type_list + + int analysis_enum = output->analysis_type_list[i]; + + /*As the domain is 2D, it is not necessary to create nodes for this analysis*/ + /*itapopo must verify if domain is not 3D. Only 2D in this version!*/ + 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,newelementslist,newx,newy,my_vertices,output->constraints); + + output->elements->Presort(); + output->nodes->Presort(); + output->vertices->Presort(); + output->loads->Presort(); + output->materials->Presort(); + output->constraints->Presort(); + + /*reset hooks for elements, loads and nodes: */ + output->elements->ResetHooks(); + output->loads->ResetHooks(); + 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; for(int i=0;inummodels;i++){ - int analysis_enum = output->analysis_type_list[i]; - //itapopo as the domain is 2D, it is not necessary to create nodes for this analysis - if(analysis_enum==StressbalanceVerticalAnalysisEnum) continue; - for(int j=0;jid=nodecounter+j+1; - newnode->sid=j; - newnode->lid=lid++; - newnode->analysis_enum=analysis_enum; - /*Initialize coord_system: Identity matrix by default*/ - for(int k=0;k<3;k++) for(int l=0;l<3;l++) newnode->coord_system[k][l]=0.0; - for(int k=0;k<3;k++) newnode->coord_system[k][k]=1.0; - /*indexing:*/ - newnode->indexingupdate=true; - Analysis* analysis=EnumToAnalysis(analysis_enum); - int *doftypes=NULL; - int numdofs=analysis->DofsPerNode(&doftypes,Domain2DhorizontalEnum,SSAApproximationEnum); - newnode->indexing.Init(numdofs,doftypes); - xDelete(doftypes); - delete analysis; - if(analysis_enum==StressbalanceAnalysisEnum) - newnode->SetApproximation(SSAApproximationEnum); - else - newnode->SetApproximation(0); - - /*Stressbalance Horiz*/ - if(analysis_enum==StressbalanceAnalysisEnum){ - // itapopo - /*Coordinate system provided, convert to coord_system matrix*/ - //XZvectorsToCoordinateSystem(&this->coord_system[0][0],&iomodel->Data(StressbalanceReferentialEnum)[j*6]); - //_assert_(sqrt( coord_system[0][0]*coord_system[0][0] + coord_system[1][0]*coord_system[1][0]) >1.e-4); - - } - output->nodes->AddObject(newnode); - } - } - - if(output->nodes->Size()) nodecounter = output->nodes->MaximumId(); - } - - /* #4: TEST THE NODES!*/ - _printf_(" Old nodes deep echo: \n"); - this->nodes->DeepEcho(); - _printf_(" New nodes deep echo: \n"); - output->nodes->DeepEcho(); - - printf(" I arrived here!!!!!!!\n"); -#ifdef _CONTINUE_NEOPZ_ - - /*Create constraints*/ - IssmDouble *spcvx=NULL; - IssmDouble *spcvy=NULL; - int numberofnodes_analysistype=this->nodes->NumberOfNodes(StressbalanceAnalysisEnum); - Vector* vspcvx=new Vector(numberofnodes_analysistype); - Vector* vspcvy=new Vector(numberofnodes_analysistype); - - IssmDouble BigNumber=1.e8; - - for(int i=0;iSetValue(i,BigNumber,INS_VAL); - vspcvy->SetValue(i,BigNumber,INS_VAL); - } - - IssmDouble absmaxspcvx=0; - IssmDouble absmaxspcvy=0; - - 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) { - vspcvx->SetValue(nodeindex,spcvalue,INS_VAL); - if(fabs(spcvalue)>absmaxspcvx) absmaxspcvx=fabs(spcvalue); - } - else { - vspcvy->SetValue(nodeindex,spcvalue,INS_VAL); - if(fabs(spcvalue)>absmaxspcvy) absmaxspcvy=fabs(spcvalue); - } - } - - /*Assemble*/ - vspcvx->Assemble(); - vspcvy->Assemble(); - - /*Serialize*/ - spcvx=vspcvx->ToMPISerial(); - spcvy=vspcvy->ToMPISerial(); - - /*Free the data*/ - delete vspcvx; - delete vspcvy; - - IssmDouble *newspcvx=NULL; - IssmDouble *newspcvy=NULL; - int *oldelements=newelementslist; //itapopo - int nods_data=numberofnodes_analysistype; - int nels_data=newnumberofelements; - int M_data=numberofnodes_analysistype; - int N_data=1; - int N_interp=newnumberofvertices;//itapopo - Options *options=NULL; - - InterpFromMeshToMesh2dx(&newspcvx,oldelements,vx,vy,nods_data,nels_data,spcvx,M_data,N_data,newx,newy,N_interp,options); - InterpFromMeshToMesh2dx(&newspcvx,oldelements,vx,vy,nods_data,nels_data,spcvx,M_data,N_data,newx,newy,N_interp,options); -// TEST: New SPCVX and VY must be tested! - - output->constraints = new Constraints(); - - nodecounter = 0; //itapopo deve começar pelo primeiro nó do StressbalanceAnalysis - int count = 0; - int constraintcounter = 0; //itapopo - IssmDouble eps = 1.e-2; - - for(int i=0;iconstraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+i+1,0,spcvx[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - } - } - count=0; - for(int i=0;iconstraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+i+1,1,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx. - count++; - } - - } - - /* #5: TEST CONSTRAINTS*/ - this->constraints->DeepEcho(); - output->constraints->DeepEcho(); - - // output->loads=static_cast(this->loads->Copy()); - // output->constraints=static_cast(this->constraints->Copy()); - // output->results=static_cast(this->results->Copy()); - - // /*reset hooks for elements, loads and nodes: */ - // output->elements->ResetHooks(); - // output->loads->ResetHooks(); - // output->materials->ResetHooks(); - - // /*do the post-processing of the datasets to get an FemModel that can actually run analyses: */ - // for(i=0;ianalysis_type_list[i]; - // output->SetCurrentConfiguration(analysis_type); - // if(i==0) VerticesDofx(output->vertices,output->parameters); //only call once, we only have one set of vertices - // SpcNodesx(output->nodes,output->constraints,output->parameters,analysis_type); - // NodesDofx(output->nodes,output->parameters,analysis_type); - // ConfigureObjectsx(output->elements,output->loads,output->nodes,output->vertices,output->materials,output->parameters); - // } - - // /*Reset current configuration: */ - // analysis_type=output->analysis_type_list[analysis_counter]; - // output->SetCurrentConfiguration(analysis_type); - - /** TODO - - - generate the required input objects for NeoPZ - AMR has fathermesh and previousmesh. On first call, these meshes are generated. - - - call NeoPZ - This creates a newmesh which will be refined. - - - get the new mesh (index,x,y) from NeoPZ - Is is doing by GetNewMesh method. - - - Create a new FemModel* that will be consistent with the new mesh - It can be done by a method. (attention with CPU #) - - - Initialize new FemModel based on new mesh (and copy the old FemModel parameters) - It can be done by a method. (attention with CPU #) - - - The last and most difficult thing to do is to update the inputs of the elements of the new mesh: - It can be done in just one method, which calls GatherInputs, InterpFromMeshToMesh and InputUpdateFromVector - - - CPU #0 will gather all inputs from FemModel - - call InterpFromMeshToMesh to interpolate them onto the new Mesh - - broadcast the new fields to all cpus - - call InputUpdateFromVector so that all inputs are updated - - - return new FemModel - - */ -#endif + analysis_type=output->analysis_type_list[i]; + output->SetCurrentConfiguration(analysis_type); + ConfigureObjectsx(output->elements,output->loads,output->nodes,output->vertices,output->materials,output->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); + } + + /*Reset current configuration: */ + analysis_type=output->analysis_type_list[this->analysis_counter]; + output->SetCurrentConfiguration(analysis_type); + + return output; +} +/*}}}*/ +void FemModel::InterpolateInputs(FemModel* newfemmodel){/*{{{*/ + + int maxinputs = MaximumNumberOfDefinitionsEnum; + + /*Figure out how many inputs we have and their respective interpolation*/ + Vector* input_interpolations=new Vector(maxinputs); + if(this->elements->Size()){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(0)); + element->GetInputsInterpolations(input_interpolations); + } + input_interpolations->Assemble(); + + /*Serialize and set output*/ + IssmDouble* input_interpolations_serial = input_interpolations->ToMPISerial(); + delete input_interpolations; + + /*Count and get enums of all inputs in old mesh*/ + int numP0inputs; + int numP1inputs; + int *P0input_enums = NULL; + int *P1input_enums = NULL; + int *P0input_interp = NULL; + int *P1input_interp = NULL; + for(int step=0;step<2;step++){ + if(step){ + P0input_enums = xNew(numP0inputs); + P1input_enums = xNew(numP1inputs); + } + numP0inputs = 0; + numP1inputs = 0; + for(int i=0;i(input_interpolations_serial[i]); + switch(inputinterp){ + case 0: + /*Input not found, go to the next*/ + break; + case P1Enum: + if(step){ + P1input_enums[numP1inputs] = i; + P1input_interp[numP1inputs] = inputinterp; + } + numP1inputs++; + break; + case P0Enum: + case DoubleInputEnum: + case IntInputEnum: + case BoolInputEnum: + if(step){ + P0input_enums[numP0inputs] = i; + P0input_interp[numP0inputs] = inputinterp; + } + numP0inputs++; + break; + default: + _error_(EnumToStringx(inputinterp)<<" Not supported yet"); + } + } + } + + printf("Found %i %i inputs\n",numP0inputs,numP1inputs); + + /*========== Deal with P0 inputs ==========*/ + int numelementsold = this->elements->NumberOfElements(); + int numelementsnew = newfemmodel->elements->NumberOfElements(); + int numverticesold = this->vertices->NumberOfVertices(); + int numverticesnew = newfemmodel->vertices->NumberOfVertices(); + IssmDouble* P0inputsold = xNew(numelementsold*numP0inputs); + IssmDouble* P0inputsnew = NULL; + + for(int i=0;i(P0inputsold); + xDelete(P0inputsnew); + + + _error_("STOP"); +} +/*}}}*/ +void FemModel::ExecuteRefinement(int &newnumberofvertices,int &newnumberofelements,IssmDouble** newx,IssmDouble** newy,IssmDouble** newz,int** pnewelementslist){/*{{{*/ + + int **newelements=NULL; + int **newsegments=NULL; + int newnumberofsegments=-1; + + /*All indexing here is in C type: 0..n-1*/ + int my_rank=IssmComm::GetRank(); + const int elementswidth=3;//just 2D mesh, tria elements + + /*Solutions which will be used to refine the elements*/ + IssmDouble* vx=NULL; //itapopo this is not being used + IssmDouble* vy=NULL; //itapopo this is not being used + IssmDouble* masklevelset=NULL; + int* newelementslist=NULL; + + this->GetMaskLevelSet(&masklevelset);//itapopo verificar se já existe um método igual a esse + + int type_process=1; //1: it refines father mesh. See AdaptiveMeshRefinement.h (.cpp) + if(my_rank==0){ + this->amr->ExecuteRefinement(type_process,vx,vy,masklevelset,newnumberofvertices,newnumberofelements,newnumberofsegments,newx,newy,newz,&newelements,&newsegments); + if(newnumberofvertices<=0 || newnumberofelements<=0 /*|| newnumberofsegments<=0*/) _error_("Error in the refinement process."); + +// printf(" PRINTING COORDINATES... \n"); + // for(int i=0;i(newnumberofelements*elementswidth); + for(int i=0;i(masklevelset); } /*}}}*/ @@ -3338,10 +3305,10 @@ /*Creating connectivity table*/ - int* connectivity=xNew(newnumberofvertices); - for(int i=0;i(newnumberofvertices); for (int i=0;i-1 && vertexidsigma=0.; newvertex->connectivity=connectivity[i]; + newvertex->clone=false;//itapopo check this vertices->AddObject(newvertex); } @@ -3373,4 +3341,5 @@ for(int i=0;iid=i+1; @@ -3387,4 +3356,5 @@ } else newtria->element_type_list=NULL; + /*Element hook*/ int matpar_id=newnumberofelements+1; //retrieve material parameter id (last pointer in femodel->materials) @@ -3392,5 +3362,5 @@ /*retrieve vertices ids*/ int* vertex_ids=xNew(elementswidth); - for(int j=0;j(newelementslist[elementswidth*i+j]); + for(int j=0;j(newelementslist[elementswidth*i+j])+1;//this Hook wants Matlab indexing /*Setting the hooks*/ newtria->numanalyses =this->nummodels; @@ -3407,4 +3377,7 @@ } } + + //itapopo there is this line in CreateElementsVerticesAndMaterials. + //elements->InputDuplicate(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum); } /*}}}*/ @@ -3424,5 +3397,224 @@ } /*}}}*/ -void FemModel::ElementsAndVerticesPartitioning(int newnumberofvertices,int newnumberofelements,int elementswidth,int* newelementslist,bool** pmy_elements,int** pmy_vertices){/*{{{*/ +void FemModel::CreateNodes(int newnumberofvertices,int* my_vertices,int nodecounter,int analysis_enum,Nodes* nodes){/*{{{*/ + + int lid=0; + + for(int j=0;jid=nodecounter+j+1; + newnode->sid=j; + newnode->lid=lid++; + newnode->analysis_enum=analysis_enum; + + /*Initialize coord_system: Identity matrix by default*/ + for(int k=0;k<3;k++) for(int l=0;l<3;l++) newnode->coord_system[k][l]=0.0; + for(int k=0;k<3;k++) newnode->coord_system[k][k]=1.0; + + /*indexing:*/ + newnode->indexingupdate=true; + + Analysis* analysis=EnumToAnalysis(analysis_enum); + int *doftypes=NULL; + int numdofs=analysis->DofsPerNode(&doftypes,Domain2DhorizontalEnum,SSAApproximationEnum); + newnode->indexing.Init(numdofs,doftypes); + xDelete(doftypes); + delete analysis; + if(analysis_enum==StressbalanceAnalysisEnum) + newnode->SetApproximation(SSAApproximationEnum); + else + newnode->SetApproximation(0); + + /*Stressbalance Horiz*/ + if(analysis_enum==StressbalanceAnalysisEnum){ + // itapopo this code is rarely used. + /*Coordinate system provided, convert to coord_system matrix*/ + //XZvectorsToCoordinateSystem(&this->coord_system[0][0],&iomodel->Data(StressbalanceReferentialEnum)[j*6]); + //_assert_(sqrt( coord_system[0][0]*coord_system[0][0] + coord_system[1][0]*coord_system[1][0]) >1.e-4); + + } + nodes->AddObject(newnode); + } + } + return; +} +/*}}}*/ +void FemModel::CreateConstraints(int newnumberofvertices,int newnumberofelements,int* newelementslist,IssmDouble* newx,IssmDouble* newy,int* my_vertices,Constraints* constraints){/*{{{*/ + + /*Get x and y of the mesh i-1*/ + int numberofvertices, numberofelements; + numberofvertices = this->vertices->NumberOfVertices(); + numberofelements = this->elements->NumberOfElements(); + + /*Get vertices coordinates*/ + IssmDouble *x = NULL; + IssmDouble *y = NULL; + IssmDouble *z = NULL; + VertexCoordinatesx(&x, &y, &z, this->vertices,false) ; + + /*Get element vertices*/ + int elementswidth = 3; //just 2D mesh in this version (just tria elements) + int* elem_vertices=xNew(elementswidth); + Vector* vid1= new Vector(numberofelements); + Vector* vid2= new Vector(numberofelements); + Vector* vid3= new Vector(numberofelements); + + /*Go through elements, and for each element, get vertices*/ + for(int i=0;ielements->Size();i++){ + Element* element=xDynamicCast(this->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 (mesh i-1)*/ + int* elementslist=NULL; + elementslist=xNew(numberofelements*elementswidth); + for(int i=0;inodes->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); + vspcvy->SetValue(i,0.,INS_VAL); + vspcvxflag->SetValue(i,0.,INS_VAL); + vspcvyflag->SetValue(i,0.,INS_VAL); + } + + 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 + vspcvx->SetValue(nodeindex,spcvalue,INS_VAL); + vspcvxflag->SetValue(nodeindex,1,INS_VAL); + } + else if(dof==1){//vy + vspcvy->SetValue(nodeindex,spcvalue,INS_VAL); + vspcvyflag->SetValue(nodeindex,1,INS_VAL); + } + else{ + /*nothing here*/ + } + } + + /*Assemble*/ + vspcvx->Assemble(); + vspcvy->Assemble(); + vspcvxflag->Assemble(); + vspcvyflag->Assemble(); + + /*Serialize*/ + spcvx=vspcvx->ToMPISerial(); + spcvy=vspcvy->ToMPISerial(); + spcvxflag=vspcvxflag->ToMPISerial(); + spcvyflag=vspcvyflag->ToMPISerial(); + /*Free the data*/ + delete vspcvx; + delete vspcvy; + delete vspcvxflag; + delete vspcvyflag; + + IssmDouble *newspcvx=NULL; + IssmDouble *newspcvy=NULL; + IssmDouble *newspcvxflag=NULL; + IssmDouble *newspcvyflag=NULL; + int nods_data=numberofvertices; + int nels_data=numberofelements; + int M_data=numberofvertices; + int N_data=1; + int N_interp=newnumberofvertices; + Options *options=new Options(); + + /*Interpolate spcvx and spcvy in the new mesh*/ + InterpFromMeshToMesh2dx(&newspcvx,elementslist,x,y,nods_data,nels_data,spcvx,M_data,N_data,newx,newy,N_interp,options); + InterpFromMeshToMesh2dx(&newspcvy,elementslist,x,y,nods_data,nels_data,spcvy,M_data,N_data,newx,newy,N_interp,options); + InterpFromMeshToMesh2dx(&newspcvxflag,elementslist,x,y,nods_data,nels_data,spcvxflag,M_data,N_data,newx,newy,N_interp,options); + InterpFromMeshToMesh2dx(&newspcvyflag,elementslist,x,y,nods_data,nels_data,spcvyflag,M_data,N_data,newx,newy,N_interp,options); + + int nodecounter = 0; //itapopo deve começar pelo primeiro nó do StressbalanceAnalysis + int count = 0; + int constraintcounter = 0; //itapopo + IssmDouble eps = 1.e-8; + + /*Now, insert the interpolated constraints in the data set (constraints)*/ + for(int i=0;i(newspcvx[i]) && newspcvxflag[i]>(1-eps)){ + constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+i+1,0,newspcvx[i],StressbalanceAnalysisEnum)); + //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + } + } + count=0; + for(int i=0;i(newspcvy[i]) && newspcvyflag[i]>(1-eps) ){ + constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+i+1,1,newspcvy[i],StressbalanceAnalysisEnum)); + //add count'th spc, on node i+1, setting dof 1 to vx. + count++; + } + + } +} +/*}}}*/ +void FemModel::UpdateElements(int newnumberofelements,int* newelementslist,bool* my_elements,int nodecounter,int analysis_counter,Elements* newelements){/*{{{*/ + + /*Update elements, set hnode. + This code is in all analysis */ + int elemcounter=0; + for(int iel=0;ielGetObjectByOffset(elemcounter); + //element update + tria->element_type_list[analysis_counter]=P1Enum; + int numnodes=3; + int* tria_node_ids=xNew(numnodes); + tria_node_ids[0]=nodecounter+newelementslist[3*iel+0]+1; //matlab indexing + tria_node_ids[1]=nodecounter+newelementslist[3*iel+1]+1; //matlab indexing + tria_node_ids[2]=nodecounter+newelementslist[3*iel+2]+1; //matlab indexing + tria->SetHookNodes(tria_node_ids,numnodes,analysis_counter); tria->nodes=NULL; + xDelete(tria_node_ids); + elemcounter++; + } + } + return; +} +/*}}}*/ +void FemModel::ElementsAndVerticesPartitioning(int& newnumberofvertices,int& newnumberofelements,int& elementswidth,int* newelementslist,bool** pmy_elements,int** pmy_vertices){/*{{{*/ int *epart=NULL; //element partitioning. @@ -3436,4 +3628,6 @@ int *my_vertices=NULL; + _assert_(newnumberofvertices>0); + _assert_(newnumberofelements>0); epart=xNew(newnumberofelements); npart=xNew(newnumberofvertices); @@ -3469,4 +3663,5 @@ will hold which vertices belong to this partition*/ for(int j=0;j