Index: /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp (revision 22240) +++ /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp (revision 22241) @@ -30,4 +30,5 @@ this->fathermesh = new TPZGeoMesh(*cp.fathermesh); this->previousmesh = new TPZGeoMesh(*cp.previousmesh); + this->refinement_type = cp.refinement_type; this->level_max = cp.level_max; this->radius_level_max = cp.radius_level_max; @@ -38,4 +39,6 @@ this->thicknesserror_threshold = cp.thicknesserror_threshold; this->deviatoricerror_threshold = cp.deviatoricerror_threshold; + this->max_deviatoricerror = cp.max_deviatoricerror; + this->max_thicknesserror = cp.max_thicknesserror; this->sid2index.clear(); this->sid2index.resize(cp.sid2index.size()); @@ -61,4 +64,5 @@ if(this->fathermesh) delete this->fathermesh; if(this->previousmesh) delete this->previousmesh; + this->refinement_type = -1; this->level_max = -1; this->radius_level_max = -1; @@ -69,4 +73,6 @@ this->thicknesserror_threshold = -1; this->deviatoricerror_threshold = -1; + this->max_deviatoricerror = -1; + this->max_thicknesserror = -1; this->sid2index.clear(); this->index2sid.clear(); @@ -79,4 +85,5 @@ this->fathermesh = NULL; this->previousmesh = NULL; + this->refinement_type = -1; this->level_max = -1; this->radius_level_max = -1; @@ -87,4 +94,6 @@ this->thicknesserror_threshold = -1; this->deviatoricerror_threshold = -1; + this->max_deviatoricerror = -1; + this->max_thicknesserror = -1; this->sid2index.clear(); this->index2sid.clear(); @@ -94,36 +103,22 @@ /*Mesh refinement methods*/ -void AdaptiveMeshRefinement::ExecuteRefinement(int numberofpoints,double* xylist,int* pnewnumberofvertices,int* pnewnumberofelements,double** px,double** py,int** pelementslist){/*{{{*/ - +void AdaptiveMeshRefinement::ExecuteRefinement(double* gl_distance,double* if_distance,double* deviatoricerror,double* thicknesserror,int* pnewnumberofvertices,int* pnewnumberofelements,double** px,double** py,int** pelementslist){/*{{{*/ + /*IMPORTANT! pelementslist are in Matlab indexing*/ /*NEOPZ works only in C indexing*/ if(!this->fathermesh || !this->previousmesh) _error_("Impossible to execute refinement: fathermesh or previousmesh is NULL!\n"); + if(this->refinement_type!=0 && this->refinement_type!=1) _error_("Impossible to execute refinement: refinement type is not defined!\n"); + + /*Input verifications*/ + if(this->deviatoricerror_threshold>0 && !deviatoricerror) _error_("deviatoricerror is NULL!\n"); + if(this->thicknesserror_threshold>0 && !thicknesserror) _error_("thicknesserror is NULL!\n"); + if(this->groundingline_distance>0 && !gl_distance) _error_("gl_distance is NULL!\n"); + if(this->icefront_distance>0 && !if_distance) _error_("if_distance is NULL!\n"); /*Intermediaries*/ bool verbose=VerboseSolution(); - - /*Refine the mesh using level max*/ - this->RefineMesh(verbose,this->previousmesh,numberofpoints,xylist); - - /*Get new geometric mesh in ISSM data structure*/ - this->GetMesh(this->previousmesh,pnewnumberofvertices,pnewnumberofelements,px,py,pelementslist); - - /*Verify the new geometry*/ - this->CheckMesh(pnewnumberofvertices,pnewnumberofelements,px,py,pelementslist); - - if(verbose) _printf_("\trefinement process done!\n"); -} -/*}}}*/ -void AdaptiveMeshRefinement::ExecuteRefinement(double* gl_elementdistance,double* if_elementdistance,double* deviatoricerror,double* thicknesserror,int* pnewnumberofvertices,int* pnewnumberofelements,double** px,double** py,int** pelementslist){/*{{{*/ - - /*IMPORTANT! pelementslist are in Matlab indexing*/ - /*NEOPZ works only in C indexing*/ - if(!this->fathermesh || !this->previousmesh) _error_("Impossible to execute refinement: fathermesh or previousmesh is NULL!\n"); - - /*Intermediaries*/ - bool verbose=true; /*Execute refinement*/ - this->RefinementProcess(verbose,gl_elementdistance,if_elementdistance,deviatoricerror,thicknesserror); + this->RefineMeshOneLevel(verbose,gl_distance,if_distance,deviatoricerror,thicknesserror); /*Get new geometric mesh in ISSM data structure*/ @@ -134,271 +129,263 @@ } /*}}}*/ -void AdaptiveMeshRefinement::RefinementProcess(bool &verbose,double* gl_elementdistance,double* if_elementdistance,double* deviatoricerror,double* thicknesserror){/*{{{*/ - - if(verbose) _printf_("\n\trefinement process started (level max = " << this->level_max << ")\n"); +void AdaptiveMeshRefinement::RefineMeshOneLevel(bool &verbose,double* gl_distance,double* if_distance,double* deviatoricerror,double* thicknesserror){/*{{{*/ /*Intermediaries*/ - TPZGeoMesh* nohangingnodesmesh=NULL; - - double mean_mask = 0; - double mean_tauerror = 0; - double mean_Herror = 0; - double group_Herror = 0; - int index,sid; - std::vector specialfatherindex; specialfatherindex.clear(); - - /*Calculate mean values*/ - /* - for(int i=0;isid2index.size();i++){ - mean_mask += masklevelset[i]; - mean_tauerror += deviatorictensorerror[i]; - mean_Herror += thicknesserror[i]; - } - mean_mask /= this->sid2index.size(); - mean_tauerror /= this->sid2index.size(); - mean_Herror /= this->sid2index.size(); - */ - if(verbose) _printf_("\t\tdeal with special elements...\n"); - /*Deal with special elements*/ - for(int i=0;ispecialelementsindex.size();i++){ - if(this->specialelementsindex[i]==-1) continue; - /*Get special element and verify*/ - TPZGeoEl* geoel=this->previousmesh->Element(this->specialelementsindex[i]); - if(!geoel)_error_("special element (sid) "<HasSubElement())_error_("special element (sid) "<NSubElements()<<" subelements!\n"); - if(geoel->MaterialId()!=this->GetElemMaterialID()) _error_("geoel->MaterialId is not GetElemMaterialID!\n"); - if(!geoel->Father())_error_("father of special element (sid) "<Father(); - TPZVec siblings; - father->GetHigherSubElements(siblings); - std::vector sidvec; sidvec.resize(siblings.size()); - for (int j=0;jindex2sid[siblings[j]->Index()]; - } - - /*Now, reset the data strucure and verify if the siblings should be deleted*/ - if(siblings.size()<4){ - /*Reset subelements in the father*/ - father->ResetSubElements(); - }else{ - if(siblings.size()!=4) _error_("element (index) "<Index()<<" has "<NSubElements()<<" subelements!\n"); - } - for (int j=0;jspecialelementsindex.size();k++){ - if(this->specialelementsindex[k]==siblings[j]->Index()){ - index = siblings[j]->Index(); - if(index<0) _error_("index is null!\n"); - sid = this->index2sid[index]; - if(sid<0) _error_("sid is null!\n"); - this->specialelementsindex[k] = -1; - this->index2sid[index] = -1; - this->sid2index[sid] = -1; - } - } - if(siblings.size()<4){ - /*Ok, the special element can be deleted*/ - siblings[j]->ResetSubElements(); - this->previousmesh->DeleteElement(siblings[j],siblings[j]->Index()); - } - } - - /*Now, verify if the father should be refined with uniform pattern (smoother)*/ - if(siblings.size()==3){//it keeps the mesh with uniform elements - /*Father has uniform subelements now*/ - TPZVec sons; - father->Divide(sons); - }else{ - specialfatherindex.push_back(father->Index()); - } - if(this->specialelementsindex[i]!=-1) _error_("special element "<previousmesh->BuildConnectivity(); - - if(verbose) _printf_("\t\tuniform refinement...\n"); - /*Deal with uniform elemnts*/ - for(int i=0;isid2index.size();i++){ - if(this->sid2index[i]==-1) continue; - /*Get element and verify*/ - TPZGeoEl* geoel=this->previousmesh->Element(this->sid2index[i]); - if(geoel->HasSubElement()) _error_("element (sid) "<NSubElements()<<" subelements!\n"); - if(geoel->MaterialId()!=this->GetElemMaterialID()) _error_("geoel->MaterialId is not GetElemMaterialID!\n"); - - /*Refine process*/ - if(thicknesserror[i]>mean_Herror) - { - int count=0; - TPZGeoEl* father=geoel->Father(); - if(father){ - for(int j=3;j<6;j++){ - index=father->Neighbour(j).Element()->Index(); - for(int k=0;k sons; - if(geoel->Level()level_max && count==0) geoel->Divide(sons); - } - else if(geoel->Level()>0) - {/*Unrefine process*/ - - /*Get siblings and verify*/ - TPZVec siblings; - geoel->Father()->GetHigherSubElements(siblings); - //if(siblings.size()<4) _error_("Impossible to refine: geoel (index) "<sid2index[i]<<" has less than 3 siblings!\n"); - if(siblings.size()>4) continue;//Father has more then 4 sons, this group should not be unrefined. - - /*Compute the error of the group*/ - group_Herror=0; - for(int j=0;jIndex(); - sid = this->index2sid[index]; - if(sid==-1) continue; - group_Herror+=thicknesserror[sid]; - } - /*Verify if this group should be unrefined*/ - if(group_Herror>0 && group_Herror<0*mean_Herror){ //itapopo - /*Reset subelements in the father*/ - this->previousmesh->Element(geoel->Father()->Index())->ResetSubElements(); - /*Delete the elements and set their indexes in the index2sid and sid2index*/ - for (int j=0;jIndex(); - sid = this->index2sid[index]; - this->index2sid[index]=-1; - if(sid!=-1) this->sid2index[sid]=-1; - this->previousmesh->DeleteElement(siblings[j],siblings[j]->Index()); - }//for j - }//if - }/*Unrefine process*/ - }//for i - this->previousmesh->BuildConnectivity(); - - if(verbose) _printf_("\t\trefine to avoid hanging nodes...\n"); - this->RefineMeshToAvoidHangingNodes(verbose,this->previousmesh); - - //nohangingnodesmesh = this->CreateRefPatternMesh(newmesh); itapopo tentar otimizar - - if(verbose) _printf_("\trefinement process done!\n"); -} -/*}}}*/ -void AdaptiveMeshRefinement::RefineMesh(bool &verbose,TPZGeoMesh* gmesh,int numberofpoints,double* xylist){/*{{{*/ - - /*Verify if there are points*/ - if(numberofpoints==0) return; - - /*Intermediaries*/ - int nelem =-1; - int side2D = 6; - double radius_h =-1; - double radius_hmax=std::max(this->radius_level_max,std::max(this->groundingline_distance,this->icefront_distance)); - int count =-1; - double mindistance=0.; - double distance =0.;; + int nelem =-1; + int side2D = 6; + int sid =-1; + int count =-1; + int criteria =-1; + int numberofcriteria =-1; + int nconformelements = this->sid2index.size(); + double gl_radius_h =-1; + double gl_radius_hmax = std::max(this->radius_level_max,this->groundingline_distance); + double if_radius_h =-1; + double if_radius_hmax = std::max(this->radius_level_max,this->icefront_distance); + double gl_groupdistance =-1; + double if_groupdistance =-1; + double d_maxerror =-1; + double t_maxerror =-1; + double deviatoric_grouperror =-1; + double thickness_grouperror =-1; + TPZGeoMesh* gmesh = NULL; TPZVec qsi(2,0.),cp(3,0.); TPZVec sons; - - /*First, delete the special elements*/ - this->DeleteSpecialElements(verbose,gmesh); + std::vector index; + + /*Calculate the number of criteria{{{*/ + numberofcriteria=0; + if(this->deviatoricerror_threshold>0) numberofcriteria++; + if(this->thicknesserror_threshold>0) numberofcriteria++; + if(this->groundingline_distance>0) numberofcriteria++; + if(this->icefront_distance>0) numberofcriteria++; + /*}}}*/ + + /*Calculate the maximum of the estimators, if requested{{{*/ + if(this->deviatoricerror_threshold>0 && this->max_deviatoricerror<0){ + for(int i=0;imax_deviatoricerror=max(this->max_deviatoricerror,deviatoricerror[i]); + } + if(this->thicknesserror_threshold>0 && this->max_thicknesserror<0){ + for(int i=0;imax_thicknesserror=max(this->max_thicknesserror,thicknesserror[i]); + } + /*}}}*/ + + /*First, verify if special elements have min distance or high errors{{{*/ + gmesh=this->previousmesh; + for(int i=0;ispecialelementsindex.size();i++){ + if(this->specialelementsindex[i]==-1) _error_("index is -1!\n"); + if(!gmesh->Element(this->specialelementsindex[i])) _error_("element is null!\n"); + if(!gmesh->Element(this->specialelementsindex[i])->Father()) _error_("father is null!\n"); + if(gmesh->Element(this->specialelementsindex[i])->HasSubElement()) _error_("special element has sub elements!\n"); + sons.clear(); + gmesh->Element(this->specialelementsindex[i])->Father()->GetHigherSubElements(sons); + /*Limits*/ + gl_radius_h = gl_radius_hmax*std::pow(this->gradation,this->level_max-gmesh->Element(this->specialelementsindex[i])->Level()); + if_radius_h = if_radius_hmax*std::pow(this->gradation,this->level_max-gmesh->Element(this->specialelementsindex[i])->Level()); + d_maxerror = this->deviatoricerror_threshold*this->max_deviatoricerror; + t_maxerror = this->thicknesserror_threshold*this->max_thicknesserror; + /*Calculate the distance and error of the group (sons)*/ + gl_groupdistance=INFINITY;if_groupdistance=INFINITY;deviatoric_grouperror=0;thickness_grouperror=0; + for(int s=0;sindex2sid[sons[s]->Index()]; + if(sid<0) continue; + if(this->groundingline_distance>0) gl_groupdistance=std::min(gl_groupdistance,gl_distance[sid]); + if(this->icefront_distance>0) if_groupdistance=std::min(if_groupdistance,if_distance[sid]); + if(this->deviatoricerror_threshold>0) deviatoric_grouperror+=deviatoricerror[sid]; + if(this->thicknesserror_threshold>0) thickness_grouperror+=thicknesserror[sid]; + } + criteria=0; + if(this->groundingline_distance>0 && gl_groupdistanceicefront_distance>0 && if_groupdistancedeviatoricerror_threshold>0 && deviatoric_grouperror>d_maxerror) criteria++; + if(this->thicknesserror_threshold>0 && thickness_grouperror>t_maxerror) criteria++; + /*Finally, it keeps the father index if it must be refine*/ + if(criteria) index.push_back(gmesh->Element(this->specialelementsindex[i])->FatherIndex()); + } + /*}}}*/ + + /*Now, detele the special elements{{{*/ + if(this->refinement_type==1) this->DeleteSpecialElements(verbose,gmesh); + else this->specialelementsindex.clear(); + /*}}}*/ + + /*Set the mesh and delete previousmesh if refinement type is 0{{{*/ + if(this->refinement_type==0){ + delete this->previousmesh; + gmesh=this->fathermesh; + } + /*}}}*/ + + /*Unrefinement process: loop over level of refinements{{{*/ + if(verbose) _printf_("\tunrefinement process...\n"); + if(verbose) _printf_("\ttotal: "); + count=0; + + nelem=gmesh->NElements();//must keep here + for(int i=0;iElement(i)) continue; + if(gmesh->Element(i)->MaterialId()!=this->GetElemMaterialID()) continue; + if(gmesh->Element(i)->HasSubElement()) continue; + if(gmesh->Element(i)->Level()==0) continue; + if(!gmesh->Element(i)->Father()) _error_("father is NULL!\n"); + /*Limits with lag*/ + gl_radius_h = this->lag*gl_radius_hmax*std::pow(this->gradation,this->level_max-gmesh->Element(i)->Level()); + if_radius_h = this->lag*if_radius_hmax*std::pow(this->gradation,this->level_max-gmesh->Element(i)->Level()); + d_maxerror = 0.05*this->max_deviatoricerror;//itapopo definir melhor + t_maxerror = 0.05*this->max_thicknesserror;//itapopo definir melhor + /*Get the sons of the father (sibilings)*/ + sons.clear(); + gmesh->Element(i)->Father()->GetHigherSubElements(sons); + if(sons.size()!=4) continue;//delete just group of 4 elements. This avoids big holes in the mesh + /*Find the minimal distance and the error of the group*/ + gl_groupdistance=INFINITY;if_groupdistance=INFINITY;deviatoric_grouperror=0;thickness_grouperror=0; + for(int s=0;sindex2sid[sons[s]->Index()]; + /*Verify if this group have solutions*/ + if(sid<0){gl_groupdistance=INFINITY;if_groupdistance=INFINITY;deviatoric_grouperror=INFINITY;thickness_grouperror=INFINITY;continue;} + /*Distance and error of the group*/ + if(this->groundingline_distance>0) gl_groupdistance=std::min(gl_groupdistance,gl_distance[sid]); + if(this->icefront_distance>0) if_groupdistance=std::min(if_groupdistance,if_distance[sid]); + if(this->deviatoricerror_threshold>0) deviatoric_grouperror+=deviatoricerror[sid]; + if(this->thicknesserror_threshold>0) thickness_grouperror+=thicknesserror[sid]; + } + /*Verify the criteria*/ + criteria=0; + if(this->groundingline_distance>0 && gl_groupdistance>gl_radius_h) criteria++; + if(this->icefront_distance>0 && if_groupdistance>if_radius_h) criteria++; + if(this->deviatoricerror_threshold>0 && deviatoric_grouperrorthicknesserror_threshold>0 && thickness_grouperrorElement(i)->Father()->ResetSubElements(); count++; + for(int s=0;sindex2sid[sons[s]->Index()]=-1;gmesh->DeleteElement(sons[s],sons[s]->Index());} + } + } + if(verbose) _printf_(""<BuildConnectivity(); + /*}}}*/ /*Refinement process: loop over level of refinements{{{*/ if(verbose) _printf_("\trefinement process (level max = "<level_max<<")\n"); - for(int h=1;h<=this->level_max;h++){ - if(verbose) _printf_("\tlevel "<NElements();//must keep here + for(int i=0;iElement(i)) continue; + if(gmesh->Element(i)->MaterialId()!=this->GetElemMaterialID()) continue; + if(gmesh->Element(i)->HasSubElement()) continue; + if(gmesh->Element(i)->Level()==this->level_max) continue; + /*Verify if this element has solutions*/ + sid=this->index2sid[gmesh->Element(i)->Index()]; + if(sid<0) continue; + /*Filter: set the region/radius for level h*/ + gl_radius_h = gl_radius_hmax*std::pow(this->gradation,this->level_max-(gmesh->Element(i)->Level()+1));//+1, current element level is gradation,this->level_max-(gmesh->Element(i)->Level()+1));//+1, current element level is deviatoricerror_threshold*this->max_deviatoricerror; + t_maxerror = this->thicknesserror_threshold*this->max_thicknesserror; + /*Verify distance and error of the element, if requested*/ + criteria=0; + if(this->groundingline_distance>0 && gl_distance[sid]icefront_distance>0 && if_distance[sid]deviatoricerror_threshold>0 && deviatoricerror[sid]>d_maxerror) criteria++; + if(this->thicknesserror_threshold>0 && thicknesserror[sid]>t_maxerror) criteria++; + /*Now, if it attends any criterion, keep the element index to refine in next step*/ + if(criteria)index.push_back(i); + } + /*Now, refine the elements*/ + for(int i=0;iElement(index[i])) DebugStop(); + if(!gmesh->Element(index[i])->HasSubElement()){gmesh->Element(index[i])->Divide(sons);count++;} + } + if(verbose) _printf_(""<BuildConnectivity(); + /*}}}*/ + + /*Now, apply smoothing and insert special elements to avoid hanging nodes{{{*/ + this->RefineMeshWithSmoothing(verbose,gmesh); + if(this->refinement_type==0) this->previousmesh=this->CreateRefPatternMesh(gmesh);//in this case, gmesh==this->fathermesh + gmesh=this->previousmesh;//previous mesh is always refined to avoid hanging nodes + this->RefineMeshToAvoidHangingNodes(verbose,gmesh); + /*}}}*/ +} +/*}}}*/ +int AdaptiveMeshRefinement::VerifyRefinementType(TPZGeoEl* geoel){/*{{{*/ + + /* + * 0 : no refinement + * 1 : special refinement (to avoid hanging nodes) + * 2 : uniform refinment + * */ + if(!geoel) _error_("geoel is NULL!\n"); + + /*Output*/ + int type=0; + int count=0; + + /*Intermediaries*/ + TPZVec sons; + + /*Loop over neighboors (sides 3, 4 and 5)*/ + for(int j=3;j<6;j++){ + sons.clear(); + geoel->Neighbour(j).Element()->GetHigherSubElements(sons); + if(sons.size()) count++; //if neighbour was refined + if(sons.size()>4) count++; //if neighbour's level is > element level+1 + } + + /*Verify and return*/ + if(count>1) type=2; + else type=count; + + return type; +} +/*}}}*/ +void AdaptiveMeshRefinement::RefineMeshWithSmoothing(bool &verbose,TPZGeoMesh* gmesh){/*{{{*/ + + /*Intermediaries*/ + int nelem =-1; + int count =-1; + int type =-1; + int typecount =-1; + + TPZVec sons; + + /*Refinement process: loop over level of refinements*/ + if(verbose) _printf_("\tsmoothing process (level max = "<level_max<<")\n"); + if(verbose) _printf_("\ttotal: "); + + count=1; + while(count>0){ count=0; - - /*Filter: set the region/radius for level h*/ - radius_h=radius_hmax*std::pow(this->gradation,this->level_max-h); - - /*Find the minimal distance of the elements (center point) to the points */ nelem=gmesh->NElements();//must keep here - for(int i=0;iElement(i)) continue; if(gmesh->Element(i)->MaterialId()!=this->GetElemMaterialID()) continue; if(gmesh->Element(i)->HasSubElement()) continue; - if(gmesh->Element(i)->Level()>=h) continue; - gmesh->Element(i)->CenterPoint(side2D,qsi); - gmesh->Element(i)->X(qsi,cp); - mindistance=radius_h; - for (int j=0;jElement(i)->Level()==this->level_max) continue; + /*loop over neighboors (sides 3, 4 and 5). Important: neighbours has the same dimension of the element*/ + type=this->VerifyRefinementType(gmesh->Element(i)); + if(type<2){ + typecount=0; + for(int j=3;j<6;j++){ + if(gmesh->Element(i)->Neighbour(j).Element()->HasSubElement()) continue; + if(gmesh->Element(i)->Neighbour(j).Element()->Index()==i) typecount++;//neighbour==this element, element at the border + if(this->VerifyRefinementType(gmesh->Element(i)->Neighbour(j).Element())==1) typecount++; + } + if(typecount>1 && type==1) type=2; + else if(typecount>2 && type==0) type=2; } - /*If the element is inside the region, refine it*/ - if(mindistanceElement(i)->Divide(sons); - count++; - } - } - if(verbose) _printf_(""<BuildConnectivity(); - /*}}}*/ - - /*Unrefinement process: loop over level of refinements{{{*/ - if(verbose) _printf_("\tunrefinement process...\n"); - for(int h=this->level_max;h>=1;h--){ - if(verbose) _printf_("\tlevel "<lag*radius_hmax*std::pow(this->gradation,this->level_max-h); - - /*Find the minimal distance of the elements (center point) to the points */ - nelem=gmesh->NElements();//must keep here - for(int i=0;iElement(i)) continue; - if(gmesh->Element(i)->MaterialId()!=this->GetElemMaterialID()) continue; - if(gmesh->Element(i)->HasSubElement()) continue; - if(gmesh->Element(i)->Level()!=h) continue; - if(!gmesh->Element(i)->Father()) _error_("father is NULL!\n"); - /*Get the sons of the father (sibilings)*/ - sons.clear(); - gmesh->Element(i)->Father()->GetHigherSubElements(sons); - if(sons.size()!=4) continue;//delete just group of 4 elements. This avoids big holes in the mesh - /*Find the minimal distance of the group*/ - mindistance=INFINITY; - for(int s=0;sCenterPoint(side2D,qsi); - sons[s]->X(qsi,cp); - for (int j=0;jradius_h){ - gmesh->Element(i)->Father()->ResetSubElements(); - count++; - for(int s=0;sDeleteElement(sons[s],sons[s]->Index()); - } - } - } - if(verbose) _printf_(""<BuildConnectivity(); - /*}}}*/ - - /*Now, insert special elements to avoid hanging nodes*/ - this->RefineMeshToAvoidHangingNodes(verbose,gmesh); -} -/*}}}*/ -void AdaptiveMeshRefinement::RefineMesh(TPZGeoMesh* gmesh,std::vector &elements){/*{{{*/ - - /*Refine elements: uniform pattern refinement*/ - for(int i=0;iElement(index); - if(geoel->HasSubElement()) _error_("Impossible to refine: geoel (index) " << index << " has subelements!\n"); - if(geoel->MaterialId()!=this->GetElemMaterialID()) _error_("Impossible to refine: geoel->MaterialId is not GetElemMaterialID!\n"); - /*Divide geoel*/ - TPZVec Sons; - geoel->Divide(Sons); - } - gmesh->BuildConnectivity(); + /*refine the element if requested*/ + if(type==2){gmesh->Element(i)->Divide(sons); count++;} + } + if(verbose){ + if(count==0) _printf_(""<BuildConnectivity(); + } } /*}}}*/ @@ -452,4 +439,5 @@ if(gmesh->Element(this->specialelementsindex[i])->Father()) gmesh->Element(this->specialelementsindex[i])->Father()->ResetSubElements(); gmesh->DeleteElement(gmesh->Element(this->specialelementsindex[i]),this->specialelementsindex[i]); + this->index2sid[this->specialelementsindex[i]]=-1; count++; } @@ -568,4 +556,5 @@ if(nvertices<=0) _error_("Impossible to create initial mesh: nvertices is <= 0!\n"); if(nelements<=0) _error_("Impossible to create initial mesh: nelements is <= 0!\n"); + if(this->refinement_type!=0 && this->refinement_type!=1) _error_("Impossible to create initial mesh: refinement type is not defined!\n"); /*Verify and creating initial mesh*/ @@ -591,23 +580,24 @@ const int mat = this->GetElemMaterialID(); TPZManVector elem(this->GetNumberOfNodes(),0); + this->index2sid.clear(); this->index2sid.resize(nelements); this->sid2index.clear(); for(int i=0;iGetNumberOfNodes();j++) elem[j]=elements[i*this->GetNumberOfNodes()+j]-1;//Convert Matlab to C indexing - /*reftype = 0: uniform, fast / reftype = 1: uniform and non-uniform (avoid hanging nodes), it is not too fast */ - const int reftype = 1; switch(this->GetNumberOfNodes()){ - case 3: this->fathermesh->CreateGeoElement(ETriangle,elem,mat,index,reftype); break; + case 3: this->fathermesh->CreateGeoElement(ETriangle,elem,mat,index,this->refinement_type); break; default: _error_("mesh not supported yet"); } /*Define the element ID*/ this->fathermesh->ElementVec()[index]->SetId(i); - /*Initialize sid2index*/ + /*Initialize sid2index and index2sid*/ this->sid2index.push_back((int)index); + this->index2sid[(int)index]=this->sid2index.size()-1;//keep the element sid } /*Build element and node connectivities*/ this->fathermesh->BuildConnectivity(); /*Set previous mesh*/ - this->previousmesh=new TPZGeoMesh(*this->fathermesh); + if(this->refinement_type==1) this->previousmesh=new TPZGeoMesh(*this->fathermesh); + else this->previousmesh=this->CreateRefPatternMesh(this->fathermesh); } /*}}}*/ @@ -622,5 +612,5 @@ int reftype = 1; long index; - + //nodes newgmesh->NodeVec().Resize(nnodes); @@ -630,9 +620,16 @@ for(int i=0;iElement(i); - TPZManVector elem(3,0); + + if(!geoel){ + index=newgmesh->ElementVec().AllocateNewElement(); + newgmesh->ElementVec()[index] = NULL; + continue; + } + + TPZManVector elem(3,0); for(int j=0;j<3;j++) elem[j] = geoel->NodeIndex(j); newgmesh->CreateGeoElement(ETriangle,elem,mat,index,reftype); - newgmesh->ElementVec()[index]->SetId(geoel->Id()); + newgmesh->ElementVec()[index]->SetId(geoel->Id()); TPZGeoElRefPattern* newgeoel = dynamic_cast*>(newgmesh->ElementVec()[index]); @@ -689,4 +686,6 @@ } } + + /*Now, build connectivities*/ newgmesh->BuildConnectivity(); @@ -715,11 +714,11 @@ //Verify if there are inf or NaN in coords for(int i=0;i<*nvertices;i++){ - if(isnan((*px)[i]) || isinf((*px)[i])) _error_("Impossible to continue: px i=" << i <<" is NaN or Inf!\n"); - if(isnan((*py)[i]) || isinf((*py)[i])) _error_("Impossible to continue: py i=" << i <<" is NaN or Inf!\n"); + if(std::isnan((*px)[i]) || std::isinf((*px)[i])) _error_("Impossible to continue: px i=" << i <<" is NaN or Inf!\n"); + if(std::isnan((*py)[i]) || std::isinf((*py)[i])) _error_("Impossible to continue: py i=" << i <<" is NaN or Inf!\n"); } for(int i=0;i<*nelements;i++){ for(int j=0;jGetNumberOfNodes();j++){ - if(isnan((*pelements)[i*GetNumberOfNodes()+j])) _error_("Impossible to continue: px i=" << i <<" is NaN!\n"); - if(isinf((*pelements)[i*GetNumberOfNodes()+j])) _error_("Impossible to continue: px i=" << i <<" is Inf!\n"); + if(std::isnan((*pelements)[i*GetNumberOfNodes()+j])) _error_("Impossible to continue: px i=" << i <<" is NaN!\n"); + if(std::isinf((*pelements)[i*GetNumberOfNodes()+j])) _error_("Impossible to continue: px i=" << i <<" is Inf!\n"); } } Index: /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h =================================================================== --- /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h (revision 22240) +++ /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h (revision 22241) @@ -54,4 +54,5 @@ * radius_h = lag * initial_radius * gradation ^ (level_max-h) */ + int refinement_type; //0 uniform (faster); 1 refpattern int level_max; //max level of refinement double radius_level_max; //initial radius which in the elements will be refined with level max @@ -63,4 +64,6 @@ double thicknesserror_threshold; //if ==0, it will not be used double deviatoricerror_threshold; //if ==0, it will not be used + double max_deviatoricerror; // Max value of the error estimator; in general, it is defined in the first time step. Attention with restart + double max_thicknesserror; // Max value of the error estimator; in general, it is defined in the first time step. Attention with restart /*}}}*/ /*Public methods{{{*/ @@ -74,5 +77,6 @@ void Initialize(); void ExecuteRefinement(int numberofpoints,double* xylist,int* pnewnumberofvertices,int* pnewnumberofelements,double** px,double** py,int** pelementslist); - void ExecuteRefinement(double* gl_elementdistance,double* if_elementdistance,double* deviatoricerror,double* thicknesserror,int* pnewnumberofvertices,int* pnewnumberofelements,double** px,double** py,int** pelementslist); + void ExecuteRefinement(int numberofpoints,double* xylist,double* deviatoricerror,double* thicknesserror,int* pnewnumberofvertices,int* pnewnumberofelements,double** px,double** py,int** pelementslist); + void ExecuteRefinement(double* gl_distance,double* if_distance,double* deviatoricerror,double* thicknesserror,int* pnewnumberofvertices,int* pnewnumberofelements,double** px,double** py,int** pelementslist); void CreateInitialMesh(int &nvertices,int &nelements,double* x,double* y,int* elements); void CheckMesh(int* nvertices,int* nelements,double** px,double** py,int** pelements); @@ -80,15 +84,14 @@ private: /*Private attributes{{{*/ - std::vector sid2index; // Vector that keeps index of PZGeoMesh elements used in the ISSM mesh (sid) - std::vector index2sid; // Vector that keeps sid of issm mesh elements used in the neopz mesh (index) - std::vector specialelementsindex; // Vector that keeps index of the special elements (created to avoid haning nodes) - TPZGeoMesh *fathermesh; // Father Mesh is the entire mesh without refinement - TPZGeoMesh *previousmesh; // Previous Mesh is the refined mesh + std::vector sid2index; // Vector that keeps index of PZGeoMesh elements used in the ISSM mesh (sid) + std::vector index2sid; // Vector that keeps sid of issm mesh elements used in the neopz mesh (index) + std::vector specialelementsindex; // Vector that keeps index of the special elements (created to avoid haning nodes) + TPZGeoMesh *fathermesh; // Entire mesh without refinement if refinement_type==1; refined with hanging nodes if efinement_type==0 + TPZGeoMesh *previousmesh; // Refined mesh without hanging nodes (it is always refpattern type), used to generate ISSM mesh /*}}}*/ /*Private methods{{{*/ - void RefinementProcess(bool &verbose,double* partiallyfloatedelements,double* masklevelset,double* deviatorictensorerror,double* thicknesserror); - void RefineMesh(bool &verbose,TPZGeoMesh* gmesh,int numberofpoints,double* xylist); - void RefineMesh(TPZGeoMesh *gmesh,std::vector &elements); - void RefineMeshToAvoidHangingNodes(bool &verbose,TPZGeoMesh* gmesh); + void RefineMeshOneLevel(bool &verbose,double* gl_distance,double* if_distance,double* deviatoricerror,double* thicknesserror); + void RefineMeshWithSmoothing(bool &verbose,TPZGeoMesh* gmesh); + void RefineMeshToAvoidHangingNodes(bool &verbose,TPZGeoMesh* gmesh); void DeleteSpecialElements(bool &verbose,TPZGeoMesh* gmesh); void GetMesh(TPZGeoMesh* gmesh,int* nvertices,int* nelements,double** px,double** py,int** pelements); @@ -98,4 +101,5 @@ void PrintGMeshVTK(TPZGeoMesh *gmesh,std::ofstream &file,bool matColor=true); int GetVTK_ElType(TPZGeoEl* gel); + int VerifyRefinementType(TPZGeoEl* geoel); /*}}}*/ }; Index: /issm/trunk-jpl/src/c/classes/AmrBamg.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/AmrBamg.cpp (revision 22240) +++ /issm/trunk-jpl/src/c/classes/AmrBamg.cpp (revision 22241) @@ -17,20 +17,21 @@ /*Constructor, copy, clean up and destructor*/ -AmrBamg::AmrBamg(IssmDouble hmin,IssmDouble hmax,int fieldenum_in,IssmDouble err_in,int keepmetric_in,IssmDouble gradation_in, - IssmDouble groundingline_resolution_in,IssmDouble groundingline_distance_in, - IssmDouble icefront_resolution_in,IssmDouble icefront_distance_in, - IssmDouble thicknesserror_resolution_in,IssmDouble thicknesserror_threshold_in, - IssmDouble deviatoricerror_resolution_in,IssmDouble deviatoricerror_threshold_in){/*{{{*/ +AmrBamg::AmrBamg(){/*{{{*/ - this->fieldenum = fieldenum_in; - this->keepmetric = keepmetric_in; - this->groundingline_resolution = groundingline_resolution_in; - this->groundingline_distance = groundingline_distance_in; - this->icefront_resolution = icefront_resolution_in; - this->icefront_distance = icefront_distance_in; - this->thicknesserror_resolution = thicknesserror_resolution_in; - this->thicknesserror_threshold = thicknesserror_threshold_in; - this->deviatoricerror_resolution = deviatoricerror_resolution_in; - this->deviatoricerror_threshold = deviatoricerror_threshold_in; + /*These attributes MUST be setup by FemModel*/ + this->fieldenum = -1;//fieldenum_in; + this->keepmetric = -1;//keepmetric_in; + this->groundingline_resolution = -1;//groundingline_resolution_in; + this->groundingline_distance = -1;//groundingline_distance_in; + this->icefront_resolution = -1;//icefront_resolution_in; + this->icefront_distance = -1;//icefront_distance_in; + this->thicknesserror_resolution = -1;//thicknesserror_resolution_in; + this->thicknesserror_threshold = -1;//thicknesserror_threshold_in; + this->thicknesserror_maximum = -1; + this->deviatoricerror_resolution = -1;//deviatoricerror_resolution_in; + this->deviatoricerror_threshold = -1;//deviatoricerror_threshold_in; + this->deviatoricerror_maximum = -1; + + /*Geometry and mesh as NULL*/ this->geometry = NULL; this->fathermesh = NULL; @@ -43,5 +44,5 @@ this->options->coeff = 1; this->options->errg = 0.1; - this->options->gradation = gradation_in; + this->options->gradation = -1; //MUST be setup by the FemModel this->options->Hessiantype = 0; this->options->maxnbv = 1e6; @@ -54,13 +55,12 @@ this->options->verbose = 0; this->options->Crack = 0; - this->options->KeepVertices = 1; /*!!!!! VERY IMPORTANT !!!!!*/ + this->options->KeepVertices = 1; /*!!!!! VERY IMPORTANT !!!!! This avoid numerical errors when remeshing*/ this->options->splitcorners = 1; - this->options->hmin = hmin; - this->options->hmax = hmax; - - this->options->err=xNew(1); - this->options->errSize[0]=1; - this->options->errSize[1]=1; - this->options->err[0] = err_in; + this->options->hmin = -1;/*MUST be setup by the FemModel*/ + this->options->hmax = -1;/*MUST be setup by the FemModel*/ + this->options->err = xNew(1); + this->options->err[0] = -1;/*MUST be setup by the FemModel*/ + this->options->errSize[0] = 1; + this->options->errSize[1] = 1; } /*}}}*/ @@ -174,2 +174,11 @@ *pelementslist = elementslist; }/*}}}*/ +void AmrBamg::SetBamgOpts(IssmDouble hmin_in,IssmDouble hmax_in,IssmDouble err_in,IssmDouble gradation_in){/*{{{*/ + + if(!this->options) _error_("AmrBamg->options is NULL!"); + + this->options->hmin = hmin_in; + this->options->hmax = hmax_in; + this->options->err[0] = err_in; + this->options->gradation= gradation_in; +}/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/AmrBamg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/AmrBamg.h (revision 22240) +++ /issm/trunk-jpl/src/c/classes/AmrBamg.h (revision 22241) @@ -21,13 +21,12 @@ IssmDouble thicknesserror_resolution; IssmDouble thicknesserror_threshold; + IssmDouble thicknesserror_maximum; IssmDouble deviatoricerror_resolution; IssmDouble deviatoricerror_threshold; + IssmDouble deviatoricerror_maximum; /* Constructor, destructor etc*/ - AmrBamg(IssmDouble hmin,IssmDouble hmax,int fieldenum_in,IssmDouble err_in,int keepmetric_in,IssmDouble gradation_in, - IssmDouble groundingline_resolution_in,IssmDouble groundingline_distance_in, - IssmDouble icefront_resolution_in,IssmDouble icefront_distance_in, - IssmDouble thicknesserror_resolution_in,IssmDouble thicknesserror_threshold_in, - IssmDouble deviatoricerror_resolution_in,IssmDouble deviatoricerror_threshold_in); + AmrBamg(); + ~AmrBamg(); @@ -35,4 +34,5 @@ void Initialize(int* elements,IssmDouble* x,IssmDouble* y,int numberofvertices,int numberofelements); void ExecuteRefinementBamg(IssmDouble* field,IssmDouble* hmaxVertices,int* pnewnumberofvertices,int *pnewnumberofelements,IssmDouble** px,IssmDouble** py,IssmDouble** pz,int** pelementslist); + void SetBamgOpts(IssmDouble hmin_in,IssmDouble hmax_in,IssmDouble err_in,IssmDouble gradation_in); private: Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 22240) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 22241) @@ -2648,6 +2648,6 @@ YC_new[i]+=Ynew[vid]; } - XC_new[i]=XC_new[i]/3; - YC_new[i]=YC_new[i]/3; + XC_new[i]=XC_new[i]/3.; + YC_new[i]=YC_new[i]/3.; } @@ -3029,5 +3029,5 @@ elementslist[elementswidth*i+0] = reCast(id1[i])+1; //InterpMesh wants Matlab indexing elementslist[elementswidth*i+1] = reCast(id2[i])+1; //InterpMesh wants Matlab indexing - elementslist[elementswidth*i+2] = reCast(id3[i])+1; //InterpMesh wants Matlab indexinf + elementslist[elementswidth*i+2] = reCast(id3[i])+1; //InterpMesh wants Matlab indexing } @@ -3608,18 +3608,21 @@ Vector* vxc = new Vector(numberofelements); Vector* vyc = new Vector(numberofelements); - IssmDouble *x = NULL; - IssmDouble *y = NULL; - IssmDouble *z = NULL; - - /*Get vertices coordinates*/ - VertexCoordinatesx(&x,&y,&z,this->vertices,false) ; + IssmDouble* x = NULL; + IssmDouble* y = NULL; + IssmDouble* z = NULL; + IssmDouble* xyz_list = NULL; + IssmDouble x1,y1,x2,y2,x3,y3; /*Insert the element center coordinates*/ for(int i=0;ielements->Size();i++){ Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); - element->GetVerticesSidList(elem_vertices); + //element->GetVerticesSidList(elem_vertices); int sid = element->Sid(); - vxc->SetValue(sid,(x[elem_vertices[0]]+x[elem_vertices[1]]+x[elem_vertices[2]])/3.,INS_VAL); - vyc->SetValue(sid,(y[elem_vertices[0]]+y[elem_vertices[1]]+y[elem_vertices[2]])/3.,INS_VAL); + element->GetVerticesCoordinates(&xyz_list); + x1 = xyz_list[3*0+0];y1 = xyz_list[3*0+1]; + x2 = xyz_list[3*1+0];y2 = xyz_list[3*1+1]; + x3 = xyz_list[3*2+0];y3 = xyz_list[3*2+1]; + vxc->SetValue(sid,(x1+x2+x3)/3.,INS_VAL); + vyc->SetValue(sid,(y1+y2+y3)/3.,INS_VAL); } @@ -3636,4 +3639,5 @@ xDelete(y); xDelete(z); + xDelete(xyz_list); xDelete(elem_vertices); delete vxc; @@ -3658,7 +3662,5 @@ int* elem_vertices = xNew(elementswidth); IssmDouble* levelset = xNew(elementswidth); - IssmDouble* x = NULL; - IssmDouble* y = NULL; - IssmDouble* z = NULL; + IssmDouble* xyz_list = NULL; Vector* vx_zerolevelset = new Vector(numberofelements); Vector* vy_zerolevelset = new Vector(numberofelements); @@ -3666,8 +3668,5 @@ IssmDouble* y_zerolevelset = NULL; int count,sid; - IssmDouble xcenter,ycenter; - - /*Get vertices coordinates*/ - VertexCoordinatesx(&x,&y,&z,this->vertices,false) ; + IssmDouble xc,yc,x1,y1,x2,y2,x3,y3; /*Use the element center coordinate if level set is zero (grounding line or ice front), otherwise set NAN*/ @@ -3676,17 +3675,22 @@ element->GetInputListOnVertices(levelset,levelset_type); element->GetVerticesSidList(elem_vertices); - sid = element->Sid(); - xcenter = NAN; - ycenter = NAN; + sid= element->Sid(); + element->GetVerticesCoordinates(&xyz_list); + x1 = xyz_list[3*0+0];y1 = xyz_list[3*0+1]; + x2 = xyz_list[3*1+0];y2 = xyz_list[3*1+1]; + x3 = xyz_list[3*2+0];y3 = xyz_list[3*2+1]; + xc = NAN; + yc = NAN; Tria* tria = xDynamicCast(element); if(tria->IsIceInElement()){/*verify if there is ice in the element*/ if(levelset[0]*levelset[1]<0. || levelset[0]*levelset[2]<0. || abs(levelset[0]*levelset[1])SetValue(sid,xcenter,INS_VAL); - vy_zerolevelset->SetValue(sid,ycenter,INS_VAL); + vx_zerolevelset->SetValue(sid,xc,INS_VAL); + vy_zerolevelset->SetValue(sid,yc,INS_VAL); + xDelete(xyz_list); } /*Assemble and serialize*/ @@ -3720,7 +3724,5 @@ xDelete(x_zerolevelset); xDelete(y_zerolevelset); - xDelete(x); - xDelete(y); - xDelete(z); + xDelete(xyz_list); delete vx_zerolevelset; delete vy_zerolevelset; @@ -4474,5 +4476,5 @@ if(this->amrbamg->deviatoricerror_threshold>0) this->GethmaxVerticesFromEstimators(hmaxvertices_serial,DeviatoricStressErrorEstimatorEnum); } - + if(my_rank==0){ this->amrbamg->ExecuteRefinementBamg(vector_serial,hmaxvertices_serial,&newnumberofvertices,&newnumberofelements,&newx,&newy,&newz,&newelementslist); @@ -4518,7 +4520,4 @@ int* elements = NULL; IssmDouble hmin,hmax,err,gradation; - IssmDouble groundingline_resolution,groundingline_distance,icefront_resolution,icefront_distance; - IssmDouble thicknesserror_resolution,thicknesserror_threshold,deviatoricerror_resolution,deviatoricerror_threshold; - int fieldenum,keepmetric; /*Get rank*/ @@ -4531,26 +4530,26 @@ this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elements); + /*Create bamg data structures for bamg*/ + this->amrbamg = new AmrBamg(); + /*Get amr parameters*/ this->parameters->FindParam(&hmin,AmrHminEnum); this->parameters->FindParam(&hmax,AmrHmaxEnum); - this->parameters->FindParam(&fieldenum,AmrFieldEnum); this->parameters->FindParam(&err,AmrErrEnum); - this->parameters->FindParam(&keepmetric,AmrKeepMetricEnum); this->parameters->FindParam(&gradation,AmrGradationEnum); - this->parameters->FindParam(&groundingline_resolution,AmrGroundingLineResolutionEnum); - this->parameters->FindParam(&groundingline_distance,AmrGroundingLineDistanceEnum); - this->parameters->FindParam(&icefront_resolution,AmrIceFrontResolutionEnum); - this->parameters->FindParam(&icefront_distance,AmrIceFrontDistanceEnum); - this->parameters->FindParam(&thicknesserror_resolution,AmrThicknessErrorResolutionEnum); - this->parameters->FindParam(&thicknesserror_threshold,AmrThicknessErrorThresholdEnum); - this->parameters->FindParam(&deviatoricerror_resolution,AmrDeviatoricErrorResolutionEnum); - this->parameters->FindParam(&deviatoricerror_threshold,AmrDeviatoricErrorThresholdEnum); - - /*Create bamg data structures for bamg*/ - this->amrbamg = new AmrBamg(hmin,hmax,fieldenum,err,keepmetric,gradation, - groundingline_resolution,groundingline_distance, - icefront_resolution,icefront_distance, - thicknesserror_resolution,thicknesserror_threshold, - deviatoricerror_resolution,deviatoricerror_threshold); + this->parameters->FindParam(&this->amrbamg->fieldenum,AmrFieldEnum); + this->parameters->FindParam(&this->amrbamg->keepmetric,AmrKeepMetricEnum); + this->parameters->FindParam(&this->amrbamg->groundingline_resolution,AmrGroundingLineResolutionEnum); + this->parameters->FindParam(&this->amrbamg->groundingline_distance,AmrGroundingLineDistanceEnum); + this->parameters->FindParam(&this->amrbamg->icefront_resolution,AmrIceFrontResolutionEnum); + this->parameters->FindParam(&this->amrbamg->icefront_distance,AmrIceFrontDistanceEnum); + this->parameters->FindParam(&this->amrbamg->thicknesserror_resolution,AmrThicknessErrorResolutionEnum); + this->parameters->FindParam(&this->amrbamg->thicknesserror_threshold,AmrThicknessErrorThresholdEnum); + this->parameters->FindParam(&this->amrbamg->thicknesserror_maximum,AmrThicknessErrorMaximumEnum); + this->parameters->FindParam(&this->amrbamg->deviatoricerror_resolution,AmrDeviatoricErrorResolutionEnum); + this->parameters->FindParam(&this->amrbamg->deviatoricerror_threshold,AmrDeviatoricErrorThresholdEnum); + this->parameters->FindParam(&this->amrbamg->deviatoricerror_maximum,AmrDeviatoricErrorMaximumEnum); + /*Set BamgOpts*/ + this->amrbamg->SetBamgOpts(hmin,hmax,err,gradation); /*Re-create original mesh and put it in bamg structure (only cpu 0)*/ @@ -4608,42 +4607,80 @@ /*Intermediaries*/ - int numberofelements = this->elements->NumberOfElements(); - IssmDouble* error_elements = NULL; - IssmDouble *x = NULL; - IssmDouble *y = NULL; - IssmDouble *z = NULL; - int *index = NULL; - IssmDouble maxerror,threshold,resolution; - int vid; - + int elementswidth = this->GetElementsWidth(); + int numberofelements = this->elements->NumberOfElements(); + int numberofvertices = this->vertices->NumberOfVertices(); + IssmDouble* maxlength = xNew(numberofelements); + IssmDouble* error_vertices = xNewZeroInit(numberofvertices); + IssmDouble* error_elements = NULL; + IssmDouble* x = NULL; + IssmDouble* y = NULL; + IssmDouble* z = NULL; + int* index = NULL; + IssmDouble maxerror,threshold,resolution,length; + IssmDouble L1,L2,L3; + int vid,v1,v2,v3; + bool refine; + + /*Fill variables*/ switch(errorestimator_type){ case ThicknessErrorEstimatorEnum: threshold = this->amrbamg->thicknesserror_threshold; resolution = this->amrbamg->thicknesserror_resolution; - this->ThicknessZZErrorEstimator(&error_elements); + maxerror = this->amrbamg->thicknesserror_maximum; + this->ThicknessZZErrorEstimator(&error_elements);//error is serial, but the calculation is parallel break; case DeviatoricStressErrorEstimatorEnum: threshold = this->amrbamg->deviatoricerror_threshold; resolution = this->amrbamg->deviatoricerror_resolution; - this->ZZErrorEstimator(&error_elements); + maxerror = this->amrbamg->deviatoricerror_maximum; + this->ZZErrorEstimator(&error_elements);//error is serial, but the calculation is parallel break; default: _error_("not implemented yet"); } - if(!error_elements) _error_("error_elements is NULL!\n"); + + /*Find the max of the estimators if it was not provided*/ + if(maxerroramrbamg->thicknesserror_maximum = maxerror;break; + case DeviatoricStressErrorEstimatorEnum: this->amrbamg->deviatoricerror_maximum = maxerror;break; + } + } /*Get mesh*/ this->GetMesh(this->vertices,this->elements,&x,&y,&z,&index); - - /*Find the max of the estimators (use error_elements)*/ - maxerror=error_elements[0]; - for(int i=0;ithreshold*maxerror){ - /*ok, fill the hmaxvertices using the element vertices*/ - for(int j=0;jGetElementsWidth();j++){ - vid=index[i*this->GetElementsWidth()+j]-1;//Matlab to C indexing + v1=index[i*elementswidth+0]-1;//Matlab to C indexing + v2=index[i*elementswidth+1]-1;//Matlab to C indexing + v3=index[i*elementswidth+2]-1;//Matlab to C indexing + L1=sqrt(pow(x[v2]-x[v1],2)+pow(y[v2]-y[v1],2)); + L2=sqrt(pow(x[v3]-x[v2],2)+pow(y[v3]-y[v2],2)); + L3=sqrt(pow(x[v1]-x[v3],2)+pow(y[v1]-y[v3],2)); + /*Fill the vectors*/ + maxlength[i] =max(L1,max(L2,L3)); + error_vertices[v1]+=error_elements[i]; + error_vertices[v2]+=error_elements[i]; + error_vertices[v3]+=error_elements[i]; + } + + /*Fill hmaxvertices with the criteria*/ + for(int i=0;i phi*maxerror*/ + if(error_elements[i]>threshold*maxerror) refine=true; + /*If the element size is closer to the resolution, verify the sum of error in the vertices*/ + if(resolution/maxlength[i]>0.85){ + for(int j=0;j0.005*maxerror) refine=true;//itapopo must be better defined + } + } + /*Now, fill the hmaxvertices if requested*/ + if(refine){ + for(int j=0;j(hmaxvertices[vid])) hmaxvertices[vid]=resolution; else hmaxvertices[vid]=min(resolution,hmaxvertices[vid]); @@ -4658,4 +4695,6 @@ xDelete(index); xDelete(error_elements); + xDelete(error_vertices); + xDelete(maxlength); } /*}}}*/ @@ -4714,7 +4753,8 @@ int my_rank = IssmComm::GetRank(); int numberofelements = this->elements->NumberOfElements(); - IssmDouble* element_label = xNewZeroInit(numberofelements); - int numberofpoints = -1; - IssmDouble* xylist = NULL; + IssmDouble* gl_distance = NULL; + IssmDouble* if_distance = NULL; + IssmDouble* deviatoricerror= NULL; + IssmDouble* thicknesserror = NULL; IssmDouble* newx = NULL; IssmDouble* newy = NULL; @@ -4724,14 +4764,14 @@ int newnumberofelements = -1; - /*Get element_label, if requested*/ - if(this->amr->groundingline_distance>0) this->GetElementLabelFromZeroLevelSet(element_label,MaskGroundediceLevelsetEnum); - if(this->amr->icefront_distance>0) this->GetElementLabelFromZeroLevelSet(element_label,MaskIceLevelsetEnum); - if(this->amr->thicknesserror_threshold>0) this->GetElementLabelFromEstimators(element_label,ThicknessErrorEstimatorEnum); - if(this->amr->deviatoricerror_threshold>0) this->GetElementLabelFromEstimators(element_label,DeviatoricStressErrorEstimatorEnum); - this->GetPointsFromElementLabel(element_label,&numberofpoints,&xylist); - + /*Get fields, if requested*/ + if(this->amr->groundingline_distance>0) this->GetElementDistanceToZeroLevelSet(&gl_distance,MaskGroundediceLevelsetEnum); + if(this->amr->icefront_distance>0) this->GetElementDistanceToZeroLevelSet(&if_distance,MaskIceLevelsetEnum); + if(this->amr->thicknesserror_threshold>0) this->ThicknessZZErrorEstimator(&thicknesserror); + if(this->amr->deviatoricerror_threshold>0) this->ZZErrorEstimator(&deviatoricerror); + if(my_rank==0){ - this->amr->ExecuteRefinement(numberofpoints,xylist,&newnumberofvertices,&newnumberofelements,&newx,&newy,&newelementslist); - newz=xNewZeroInit(newnumberofvertices); + this->amr->ExecuteRefinement(gl_distance,if_distance,deviatoricerror,thicknesserror, + &newnumberofvertices,&newnumberofelements,&newx,&newy,&newelementslist); + newz=xNewZeroInit(newnumberofvertices); if(newnumberofvertices<=0 || newnumberofelements<=0) _error_("Error in the ReMeshNeopz."); } @@ -4760,6 +4800,8 @@ /*Cleanup*/ - xDelete(element_label); - xDelete(xylist); + xDelete(deviatoricerror); + xDelete(thicknesserror); + xDelete(gl_distance); + xDelete(if_distance); } /*}}}*/ @@ -4805,4 +4847,5 @@ this->SetRefPatterns(); this->amr = new AdaptiveMeshRefinement(); + this->amr->refinement_type = 1;//1 is refpattern; 0 is uniform (faster) this->amr->level_max = level_max; this->amr->radius_level_max = radius_level_max; @@ -4824,138 +4867,4 @@ } /*}}}*/ -void FemModel::GetPointsFromElementLabel(IssmDouble* element_label,int* pnumberofpoints,IssmDouble** pxylist){/*{{{*/ - - if(!element_label) _error_("element_label is NULL!\n"); - - /*Outputs*/ - int numberofpoints = -1; - IssmDouble* xylist = NULL; - - /*Intermediaries*/ - int numberofelements = this->elements->NumberOfElements(); - int count = -1; - IssmDouble* xc = NULL; - IssmDouble* yc = NULL; - - /*First, find the number of labeled elements (points)*/ - count=0; - for(int i=0;iDBL_EPSILON) count++; - } - - /*Set number of points*/ - numberofpoints=count; - if(count>0) xylist=xNew(2*numberofpoints); - - /*Get element center coordinates*/ - this->GetElementCenterCoordinates(&xc,&yc); - - /*Now, fill xylist data*/ - count=0; - for(int i=0;iDBL_EPSILON){ - xylist[2*count] = xc[i]; - xylist[2*count+1] = yc[i]; - count++; - } - } - - /*Assign pointers*/ - (*pxylist)=xylist; - *pnumberofpoints=numberofpoints; - - /*Cleanup*/ - xDelete(xc); - xDelete(yc); -} -/*}}}*/ -void FemModel::GetElementLabelFromZeroLevelSet(IssmDouble* element_label,int levelset_type){/*{{{*/ - - /*Here, "zero level set" means grounding line or ice front, depending on the level set type*/ - /*element_label is 1 if the element zero level set, NAN otherwise*/ - if(levelset_type!=MaskGroundediceLevelsetEnum && levelset_type!=MaskIceLevelsetEnum) _error_("level set type not implemented yet!"); - if(!element_label) _error_("element_label is NULL!\n"); - - /*Intermediaries*/ - int elementswidth = this->GetElementsWidth(); - int numberofelements = this->elements->NumberOfElements(); - int* elem_vertices = xNew(elementswidth); - IssmDouble* levelset = xNew(elementswidth); - Vector* velement_label = new Vector(numberofelements); - IssmDouble* element_label_serial = NULL; - int sid = -1; - IssmDouble label = -1.; - - /*Use the element center coordinate if level set is zero (grounding line or ice front), otherwise set NAN*/ - for(int i=0;ielements->Size();i++){ - Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); - element->GetInputListOnVertices(levelset,levelset_type); - element->GetVerticesSidList(elem_vertices); - sid = element->Sid(); - label = NAN; - Tria* tria = xDynamicCast(element); - if(tria->IsIceInElement()){/*verify if there is ice in the element*/ - if(levelset[0]*levelset[1]<0. || levelset[0]*levelset[2]<0. || - abs(levelset[0]*levelset[1])SetValue(sid,label,INS_VAL); - } - - /*Assemble and serialize*/ - velement_label->Assemble(); - element_label_serial=velement_label->ToMPISerial(); - - /*Merge with the output*/ - for(int i=0;i(element_label_serial[i])) element_label[i]=element_label_serial[i]; - else; //do nothing - } - - /*Cleanup*/ - xDelete(elem_vertices); - xDelete(levelset); - xDelete(element_label_serial); - delete velement_label; -} -/*}}}*/ -void FemModel::GetElementLabelFromEstimators(IssmDouble* element_label,int estimator_type){/*{{{*/ - - /*element_label is 1 if the element zero level set, NAN otherwise*/ - if(!element_label) _error_("element_label is NULL!\n"); - - /*Intermediaries*/ - int numberofelements = this->elements->NumberOfElements(); - IssmDouble* elementerror = NULL; - IssmDouble threshold = -1.; - IssmDouble maxerror = -1.; - - switch(estimator_type){ - case ThicknessErrorEstimatorEnum: - threshold=this->amr->thicknesserror_threshold; - this->ThicknessZZErrorEstimator(&elementerror); - break; - case DeviatoricStressErrorEstimatorEnum: - threshold=this->amr->deviatoricerror_threshold; - this->ZZErrorEstimator(&elementerror); - break; - default: _error_("not implemented yet"); - } - - /*Find the max of the estimators*/ - maxerror=elementerror[0]; - for(int i=0;ithreshold*maxerror) element_label[i]=1.; - else; //do nothing - } - - /*Cleanup*/ - xDelete(elementerror); -} -/*}}}*/ void FemModel::GetElementDistanceToZeroLevelSet(IssmDouble** pelementdistance,int levelset_type){/*{{{*/ @@ -4970,34 +4879,44 @@ /*Intermediaries*/ - int numberofelements = this->elements->NumberOfElements(); - IssmDouble* levelset_points= NULL; - IssmDouble* xc = NULL; - IssmDouble* yc = NULL; + int numberofelements = this->elements->NumberOfElements(); + Vector* velementdistance = new Vector(numberofelements); + IssmDouble* levelset_points = NULL; + IssmDouble* xyz_list = NULL; + IssmDouble mindistance,distance; + IssmDouble xc,yc,x1,y1,x2,y2,x3,y3; int numberofpoints; - IssmDouble distance; - - /*Get element center coordinates*/ - this->GetElementCenterCoordinates(&xc,&yc); - - /*Get points which level set is zero (center of elements with zero level set)*/ + + /*Get points which level set is zero (center of elements with zero level set, levelset_points is serial)*/ this->GetZeroLevelSetPoints(&levelset_points,numberofpoints,levelset_type); - /*Find the minimal element distance to the zero levelset (grounding line or ice front)*/ - elementdistance=xNew(numberofelements); - for(int i=0;ielements->Size();i++){//parallel + Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); + int sid = element->Sid(); + element->GetVerticesCoordinates(&xyz_list); + x1 = xyz_list[3*0+0];y1 = xyz_list[3*0+1]; + x2 = xyz_list[3*1+0];y2 = xyz_list[3*1+1]; + x3 = xyz_list[3*2+0];y3 = xyz_list[3*2+1]; + xc = (x1+x2+x3)/3.; + yc = (y1+y2+y3)/3.; + mindistance=INFINITY; + /*Loop over each point (where level set is zero)*/ for(int j=0;jSetValue(sid,mindistance,INS_VAL); + xDelete(xyz_list); } + /*Assemble*/ + velementdistance->Assemble(); + /*Assign the pointer*/ - (*pelementdistance)=elementdistance; + (*pelementdistance)=velementdistance->ToMPISerial(); /*Cleanup*/ xDelete(levelset_points); - xDelete(xc); - xDelete(yc); + xDelete(xyz_list); + delete velementdistance; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 22240) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 22241) @@ -193,7 +193,4 @@ void ReMeshNeopz(int* pnewnumberofvertices,int* pnewnumberofelements,IssmDouble** pnewx,IssmDouble** pnewy,IssmDouble** pnewz,int** pnewelementslist); void InitializeAdaptiveRefinementNeopz(void); - void GetPointsFromElementLabel(IssmDouble* element_label,int* numberofpoints,IssmDouble** xylist); - void GetElementLabelFromZeroLevelSet(IssmDouble* element_label,int levelset_type); - void GetElementLabelFromEstimators(IssmDouble* element_label,int estimator_type); void GetElementDistanceToZeroLevelSet(IssmDouble** pelementdistance,int levelset_type); void SetRefPatterns(void); Index: /issm/trunk-jpl/src/c/cores/transient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 22240) +++ /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 22241) @@ -26,5 +26,5 @@ int output_frequency; int recording_frequency; - int domaintype,groundingline_migration,smb_model,amr_frequency; + int domaintype,groundingline_migration,smb_model,amr_frequency,amr_restart; int numoutputs; Analysis *analysis = NULL; @@ -65,7 +65,9 @@ if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,TransientRequestedOutputsEnum); - #ifdef _HAVE_NEOPZ_ - bool ismismip = false;//itapopo testing restart - if(ismismip) femmodel->ReMesh(); + #ifdef _HAVE_BAMG_ //#ifdef _HAVE_NEOPZ_ itapopo + if(amr_frequency){ + femmodel->parameters->FindParam(&amr_restart,AmrRestartEnum); + if(amr_restart) femmodel->ReMesh(); + } #endif Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 22240) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 22241) @@ -136,5 +136,8 @@ parameters->AddObject(iomodel->CopyConstantObject("md.amr.icefront_distance",AmrIceFrontDistanceEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.amr.thicknesserror_threshold",AmrThicknessErrorThresholdEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.thicknesserror_maximum",AmrThicknessErrorMaximumEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.amr.deviatoricerror_threshold",AmrDeviatoricErrorThresholdEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.deviatoricerror_maximum",AmrDeviatoricErrorMaximumEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.restart",AmrRestartEnum)); break; #endif @@ -153,6 +156,9 @@ parameters->AddObject(iomodel->CopyConstantObject("md.amr.thicknesserror_resolution",AmrThicknessErrorResolutionEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.amr.thicknesserror_threshold",AmrThicknessErrorThresholdEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.thicknesserror_maximum",AmrThicknessErrorMaximumEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.amr.deviatoricerror_resolution",AmrDeviatoricErrorResolutionEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.amr.deviatoricerror_threshold",AmrDeviatoricErrorThresholdEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.deviatoricerror_maximum",AmrDeviatoricErrorMaximumEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.restart",AmrRestartEnum)); /*Convert fieldname to enum and put it in params*/ iomodel->FindConstant(&fieldname,"md.amr.fieldname"); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 22240) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 22241) @@ -870,4 +870,5 @@ TransientAmrFrequencyEnum, AmrTypeEnum, + AmrRestartEnum, AmrNeopzEnum, AmrLevelMaxEnum, @@ -887,6 +888,8 @@ AmrThicknessErrorResolutionEnum, AmrThicknessErrorThresholdEnum, + AmrThicknessErrorMaximumEnum, AmrDeviatoricErrorResolutionEnum, AmrDeviatoricErrorThresholdEnum, + AmrDeviatoricErrorMaximumEnum, DeviatoricStressErrorEstimatorEnum, ThicknessErrorEstimatorEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 22240) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 22241) @@ -843,4 +843,5 @@ case TransientAmrFrequencyEnum : return "TransientAmrFrequency"; case AmrTypeEnum : return "AmrType"; + case AmrRestartEnum : return "AmrRestart"; case AmrNeopzEnum : return "AmrNeopz"; case AmrLevelMaxEnum : return "AmrLevelMax"; @@ -860,6 +861,8 @@ case AmrThicknessErrorResolutionEnum : return "AmrThicknessErrorResolution"; case AmrThicknessErrorThresholdEnum : return "AmrThicknessErrorThreshold"; + case AmrThicknessErrorMaximumEnum : return "AmrThicknessErrorMaximum"; case AmrDeviatoricErrorResolutionEnum : return "AmrDeviatoricErrorResolution"; case AmrDeviatoricErrorThresholdEnum : return "AmrDeviatoricErrorThreshold"; + case AmrDeviatoricErrorMaximumEnum : return "AmrDeviatoricErrorMaximum"; case DeviatoricStressErrorEstimatorEnum : return "DeviatoricStressErrorEstimator"; case ThicknessErrorEstimatorEnum : return "ThicknessErrorEstimator"; Index: /issm/trunk-jpl/src/m/classes/amr.js =================================================================== --- /issm/trunk-jpl/src/m/classes/amr.js (revision 22240) +++ /issm/trunk-jpl/src/m/classes/amr.js (revision 22241) @@ -19,6 +19,8 @@ this.thicknesserror_resolution = 500; this.thicknesserror_threshold = 0; + this.thicknesserror_maximum = 0; this.deviatoricerror_resolution = 500; this.deviatoricerror_threshold = 0; + this.deviatoricerror_maximum = 0; }// }}} this.disp= function(){// {{{ @@ -35,6 +37,8 @@ fielddisplay(this,'thicknesserror_resolution','element length when thickness error estimator is used'); fielddisplay(this,'thicknesserror_threshold','maximum threshold thickness error permitted'); + fielddisplay(this,'thicknesserror_maximum','maximum thickness error permitted'); fielddisplay(this,'deviatoricerror_resolution','element length when deviatoric stress error estimator is used'); fielddisplay(this,'deviatoricerror_threshold','maximum threshold deviatoricstress error permitted'); + fielddisplay(this,'deviatoricerror_maximum','maximum deviatoricstress error permitted'); }// }}} this.classname= function(){// {{{ @@ -53,6 +57,8 @@ checkfield(md,'fieldname','amr.thicknesserror_resolution','numel',[1],'>',0,'<',this.hmax,'NaN',1); checkfield(md,'fieldname','amr.thicknesserror_threshold','numel',[1],'>=',0,'<=',1,'NaN',1); + checkfield(md,'fieldname','amr.thicknesserror_maximum','numel',[1],'>=',0,'NaN',1,'Inf',1); checkfield(md,'fieldname','amr.deviatoricerror_resolution','numel',[1],'>',0,'<',this.hmax,'NaN',1); checkfield(md,'fieldname','amr.deviatoricerror_threshold','numel',[1],'>=',0,'<=',1,'NaN',1); + checkfield(md,'fieldname','amr.deviatoricerror_maximum','numel',[1],'>=',0,'NaN',1,'Inf',1); } // }}} this.marshall=function(md,prefix,fid) { //{{{ @@ -70,6 +76,8 @@ WriteData(fid,prefix,'object',this,'fieldname','thicknesserror_resolution','format','Double'); WriteData(fid,prefix,'object',this,'fieldname','thicknesserror_threshold','format','Double'); + WriteData(fid,prefix,'object',this,'fieldname','thicknesserror_maximum','format','Double'); WriteData(fid,prefix,'object',this,'fieldname','deviatoricerror_resolution','format','Double'); WriteData(fid,prefix,'object',this,'fieldname','deviatoricerror_threshold','format','Double'); + WriteData(fid,prefix,'object',this,'fieldname','deviatoricerror_maximum','format','Double'); }//}}} this.fix=function() { //{{{ @@ -89,6 +97,8 @@ this.thicknesserror_resolution = 0.; this.thicknesserror_threshold = 0.; + this.thicknesserror_maximum = 0.; this.deviatoricerror_resolution = 0.; this.deviatoricerror_threshold = 0.; + this.deviatoricerror_maximum = 0.; this.setdefaultparameters(); Index: /issm/trunk-jpl/src/m/classes/amr.m =================================================================== --- /issm/trunk-jpl/src/m/classes/amr.m (revision 22240) +++ /issm/trunk-jpl/src/m/classes/amr.m (revision 22241) @@ -18,6 +18,9 @@ thicknesserror_resolution = 0.; thicknesserror_threshold = 0.; + thicknesserror_maximum = 0.; deviatoricerror_resolution = 0.; deviatoricerror_threshold = 0.; + deviatoricerror_maximum = 0.; + restart=0.; end methods (Static) @@ -86,6 +89,11 @@ self.thicknesserror_resolution=500.; self.thicknesserror_threshold=0.; + self.thicknesserror_maximum=0.; self.deviatoricerror_resolution=500.; self.deviatoricerror_threshold=0.; + self.deviatoricerror_maximum=0.; + + %is restart? This calls femmodel->ReMesh() before first time step. + self.restart=0; end % }}} @@ -103,6 +111,9 @@ md = checkfield(md,'fieldname','amr.thicknesserror_resolution','numel',[1],'>',0,'<',self.hmax,'NaN',1); md = checkfield(md,'fieldname','amr.thicknesserror_threshold','numel',[1],'>=',0,'<=',1,'NaN',1); + md = checkfield(md,'fieldname','amr.thicknesserror_maximum','numel',[1],'>=',0,'NaN',1,'Inf',1); md = checkfield(md,'fieldname','amr.deviatoricerror_resolution','numel',[1],'>',0,'<',self.hmax,'NaN',1); md = checkfield(md,'fieldname','amr.deviatoricerror_threshold','numel',[1],'>=',0,'<=',1,'NaN',1); + md = checkfield(md,'fieldname','amr.deviatoricerror_maximum','numel',[1],'>=',0,'NaN',1,'Inf',1); + md = checkfield(md,'fieldname','amr.restart','numel',[1],'>=',0,'<=',1,'NaN',1); end % }}} function disp(self) % {{{ @@ -120,7 +131,9 @@ fielddisplay(self,'thicknesserror_resolution',['element length when thickness error estimator is used']); fielddisplay(self,'thicknesserror_threshold',['maximum threshold thickness error permitted']); + fielddisplay(self,'thicknesserror_maximum',['maximum thickness error permitted']); fielddisplay(self,'deviatoricerror_resolution',['element length when deviatoric stress error estimator is used']); fielddisplay(self,'deviatoricerror_threshold',['maximum threshold deviatoricstress error permitted']); - + fielddisplay(self,'deviatoricerror_maximum',['maximum deviatoricstress error permitted']); + fielddisplay(self,'restart',['indicates if ReMesh() will call before first time step']); end % }}} function marshall(self,prefix,md,fid) % {{{ @@ -139,6 +152,9 @@ WriteData(fid,prefix,'object',self,'class','amr','fieldname','thicknesserror_resolution','format','Double'); WriteData(fid,prefix,'object',self,'class','amr','fieldname','thicknesserror_threshold','format','Double'); + WriteData(fid,prefix,'object',self,'class','amr','fieldname','thicknesserror_maximum','format','Double'); WriteData(fid,prefix,'object',self,'class','amr','fieldname','deviatoricerror_resolution','format','Double'); WriteData(fid,prefix,'object',self,'class','amr','fieldname','deviatoricerror_threshold','format','Double'); + WriteData(fid,prefix,'object',self,'class','amr','fieldname','deviatoricerror_maximum','format','Double'); + WriteData(fid,prefix,'object',self,'class','amr','fieldname','restart','format','Integer'); end % }}} Index: /issm/trunk-jpl/src/m/classes/amr.py =================================================================== --- /issm/trunk-jpl/src/m/classes/amr.py (revision 22240) +++ /issm/trunk-jpl/src/m/classes/amr.py (revision 22241) @@ -24,6 +24,8 @@ self.thicknesserror_resolution = 0. self.thicknesserror_threshold = 0. + self.thicknesserror_maximum = 0. self.deviatoricerror_resolution= 0. self.deviatoricerror_threshold = 0. + self.deviatoricerror_maximum = 0. #set defaults self.setdefaultparameters() @@ -42,6 +44,8 @@ string="%s\n%s"%(string,fielddisplay(self,"thicknesserror_resolution","element length when thickness error estimator is used")) string="%s\n%s"%(string,fielddisplay(self,"thicknesserror_threshold","maximum threshold thickness error permitted")) + string="%s\n%s"%(string,fielddisplay(self,"thicknesserror_maximum","maximum thickness error permitted")) string="%s\n%s"%(string,fielddisplay(self,"deviatoricerror_resolution","element length when deviatoric stress error estimator is used")) string="%s\n%s"%(string,fielddisplay(self,"deviatoricerror_threshold","maximum threshold deviatoricstress error permitted")) + string="%s\n%s"%(string,fielddisplay(self,"deviatoricerror_maximum","maximum deviatoricstress error permitted")) return string #}}} @@ -59,6 +63,8 @@ self.thicknesserror_resolution = 500. self.thicknesserror_threshold = 0 + self.thicknesserror_maximum = 0 self.deviatoricerror_resolution= 500. self.deviatoricerror_threshold = 0 + self.deviatoricerror_maximum = 0 return self #}}} @@ -74,6 +80,8 @@ md = checkfield(md,'fieldname','amr.thicknesserror_resolution','numel',[1],'>',0,'<',self.hmax,'NaN',1); md = checkfield(md,'fieldname','amr.thicknesserror_threshold','numel',[1],'>=',0,'<=',1,'NaN',1); + md = checkfield(md,'fieldname','amr.thicknesserror_maximum','numel',[1],'>=',0,'NaN',1,'Inf',1); md = checkfield(md,'fieldname','amr.deviatoricerror_resolution','numel',[1],'>',0,'<',self.hmax,'NaN',1); md = checkfield(md,'fieldname','amr.deviatoricerror_threshold','numel',[1],'>=',0,'<=',1,'NaN',1); + md = checkfield(md,'fieldname','amr.deviatoricerror_maximum','numel',[1],'>=',0,'NaN',1,'Inf',1); return md # }}} @@ -92,5 +100,7 @@ WriteData(fid,prefix,'object',self,'fieldname','thicknesserror_resolution','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','thicknesserror_threshold','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','thicknesserror_maximum','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','deviatoricerror_resolution','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','deviatoricerror_threshold','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','deviatoricerror_maximum','format','Double'); # }}} Index: /issm/trunk-jpl/src/m/contrib/tsantos/AMRexportVTK.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/tsantos/AMRexportVTK.m (revision 22240) +++ /issm/trunk-jpl/src/m/contrib/tsantos/AMRexportVTK.m (revision 22241) @@ -60,8 +60,17 @@ timestep=step; - - points=[model.results.TransientSolution(step).MeshX model.results.TransientSolution(step).MeshY zeros(size(model.results.TransientSolution(step).MeshX))]; + if(isfield(model.results.TransientSolution,'MeshElements')) + index = model.results.TransientSolution(step).MeshElements; + x = model.results.TransientSolution(step).MeshX; + y = model.results.TransientSolution(step).MeshY; + else + index = model.mesh.elements; + x = model.mesh.x; + y = model.mesh.y; + end + + points=[x y zeros(size(x))]; [num_of_points,dim]=size(points); - [num_of_elt]=size(model.results.TransientSolution(step).MeshElements,1); + [num_of_elt]=size(index,1); fid = fopen(strcat(path,filesep,name,filesep,'timestep.vtk',int2str(timestep),'.vtk'),'w+'); @@ -86,5 +95,5 @@ end s=cell2mat(horzcat(s,{'\n'})); - fprintf(fid,s,[(point_per_elt)*ones(num_of_elt,1) model.results.TransientSolution(step).MeshElements-1]'); + fprintf(fid,s,[(point_per_elt)*ones(num_of_elt,1) index-1]'); fprintf(fid,'CELL_TYPES %d\n',num_of_elt); Index: /issm/trunk-jpl/src/m/contrib/tsantos/mismip/gl_position.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/tsantos/mismip/gl_position.m (revision 22240) +++ /issm/trunk-jpl/src/m/contrib/tsantos/mismip/gl_position.m (revision 22241) @@ -3,7 +3,13 @@ %initialization of some variables data = md.results.TransientSolution(step).MaskGroundediceLevelset; - index = md.results.TransientSolution(step).MeshElements; - x = md.results.TransientSolution(step).MeshX; - y = md.results.TransientSolution(step).MeshY; + if(isfield(md.results.TransientSolution,'MeshElements')) + index = md.results.TransientSolution(step).MeshElements; + x = md.results.TransientSolution(step).MeshX; + y = md.results.TransientSolution(step).MeshY; + else + index = md.mesh.elements; + x = md.mesh.x; + y = md.mesh.y; + end numberofelements = size(index,1); elementslist = 1:numberofelements; Index: /issm/trunk-jpl/src/m/contrib/tsantos/mismip/ice_evolution.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/tsantos/mismip/ice_evolution.m (revision 22240) +++ /issm/trunk-jpl/src/m/contrib/tsantos/mismip/ice_evolution.m (revision 22241) @@ -2,32 +2,90 @@ %iv: ice volume %ivaf: ice volume above floatation -%GLy40 : grounding line position @ y=40km +%GL_y : grounding line position @ y (y comes in m) %nelem : number of elements +% usage: +% +% Default: y=40km, i0=1 +% [ga iv ivaf GL_y nelem t] = ice_evolution(md); +% +% Default: y=40km +% [ga iv ivaf GL_y nelem t] = ice_evolution(md,i0); +% +% Use this for y the borders (y=0 or y=ymax) +% [ga iv ivaf GL_y nelem t] = ice_evolution(md,i0,y); +% +% -function [ga iv ivaf GLy40 nelem t] = ice_evolution(md), +function [ga iv ivaf GL_y nelem t] = ice_evolution(varargin), ga = []; iv = []; ivaf = []; - GLy40 = []; + GL_y = []; nelem = []; t = []; + + if(nargin==0) + error('it is necessary the model!') + elseif(nargin==1) + % Defoult is y=40km + i0 = 1; + y0 = 35000; + y1 = 45000; + dy = 100; + y = 40000; + elseif(nargin==2) + i0=varargin{2}; + % Defoult is y=40km + y0 = 35000; + y1 = 45000; + dy = 100; + y = 40000; + elseif(nargin==3) + i0 = varargin{2}; + y = varargin{3}; + dy = 10; + y0 = y-dy; + y1 = y+dy; + else + error('number of inputs is more than 3'); + end + %set the model + md =varargin{1}; nsteps = length(md.results.TransientSolution); - for i=1:nsteps, + + for i=i0:nsteps, ga(i) = md.results.TransientSolution(i).GroundedArea; iv(i) = md.results.TransientSolution(i).IceVolume; ivaf(i) = md.results.TransientSolution(i).IceVolumeAboveFloatation; - nelem(i) = size(md.results.TransientSolution(i).MeshElements,1); + if(isfield(md.results.TransientSolution,'MeshElements')) + nelem(i) = size(md.results.TransientSolution(i).MeshElements,1); + else + nelem(i) = md.mesh.numberofelements; + end t(i) = md.results.TransientSolution(i).time; - %find GL position at y=40km + %find GL position between y0 and y1 [glx gly] = gl_position(md,i,0); - pos=find(gly<45000 & gly > 35000); - x=gly(pos); - v=glx(pos); - xq=[38000:100:42000]; - vq = interp1(x,v,xq,'linear'); - pos=find(xq==40000); - GLy40(i)=vq(pos); + pos = find(glyy0); + x = gly(pos); + v = glx(pos); + if(length(pos)==0) + error('pos is null') + elseif(length(pos)==1) + %this should be used for y=0 or y=ymax + GL_y(i) = v; + else + %this should be used when y is inside the domain; so, use linear interpolation + xq = [y0:dy:y1]; + vq = interp1(x,v,xq,'linear'); + pos = find(xq==y); + if(pos) + GL_y(i) = vq(pos); + else + error('pos is null') + end + end + end Index: /issm/trunk-jpl/src/m/contrib/tsantos/remesh.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/tsantos/remesh.m (revision 22240) +++ /issm/trunk-jpl/src/m/contrib/tsantos/remesh.m (revision 22241) @@ -58,5 +58,5 @@ NewModel.geometry.surface = md.results.TransientSolution(end).Surface; NewModel.geometry.base = md.results.TransientSolution(end).Base; -%NewModel.geometry.bed = md.geometry.bed; %use from parameterize +NewModel.geometry.bed = md.results.TransientSolution(end).Bed;%md.geometry.bed; %use from parameterize NewModel.geometry.thickness = md.results.TransientSolution(end).Thickness; NewModel.mask.groundedice_levelset = md.results.TransientSolution(end).MaskGroundediceLevelset; @@ -72,4 +72,5 @@ NewModel.cluster = md.cluster; NewModel.transient = md.transient; +NewModel.amr = md.amr; mdOut = NewModel;