Index: /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp (revision 21918) +++ /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp (revision 21919) @@ -10,8 +10,5 @@ #include "./AdaptiveMeshRefinement.h" -#include "TPZVTKGeoMesh.h" -#include "../shared/shared.h" -#include "pzgeotriangle.h" -#include "pzreftriangle.h" + using namespace pzgeom; @@ -40,23 +37,15 @@ this->sid2index.resize(cp.sid2index.size()); for(int i=0;isid2index[i]=cp.sid2index[i]; - + this->index2sid.clear(); + this->index2sid.resize(cp.index2sid.size()); + for(int i=0;iindex2sid[i]=cp.index2sid[i]; + this->specialelementsindex.clear(); + this->specialelementsindex.resize(cp.specialelementsindex.size()); + for(int i=0;ispecialelementsindex[i]=cp.specialelementsindex[i]; + return *this; - } /*}}}*/ AdaptiveMeshRefinement::~AdaptiveMeshRefinement(){/*{{{*/ - - bool ismismip = false; - if(ismismip){//itapopo - TPZFileStream fstr; - std::stringstream ss; - - ss << this->levelmax; - std::string AMRfile = "/home/santos/L" + ss.str() + "_amr.txt"; - - fstr.OpenWrite(AMRfile.c_str()); - int withclassid = 1; - this->Write(fstr,withclassid); - } this->CleanUp(); gRefDBase.clear(); @@ -73,4 +62,6 @@ this->regionlevelmax = -1; this->sid2index.clear(); + this->index2sid.clear(); + this->specialelementsindex.clear(); } /*}}}*/ @@ -84,79 +75,9 @@ this->regionlevel1 = -1; this->regionlevelmax = -1; + this->step = 0; + this->smooth_frequency = 1; this->sid2index.clear(); -} -/*}}}*/ -int AdaptiveMeshRefinement::ClassId() const{/*{{{*/ - return 13829430; //Antartic area with ice shelves (km^2) -} -/*}}}*/ -void AdaptiveMeshRefinement::Read(TPZStream &buf,void *context){/*{{{*/ - - try - { - /* Read the id context*/ - TPZSaveable::Read(buf,context); - /* Read class id*/ - int classid; - buf.Read(&classid,1); - /* Verify the class id*/ - if(classid!=this->ClassId()) _error_("AdaptiveMeshRefinement::Read: Error in restoring AdaptiveMeshRefinement!\n"); - /* Read simple attributes */ - buf.Read(&this->levelmax,1); - buf.Read(&this->elementswidth,1); - buf.Read(&this->regionlevel1,1); - buf.Read(&this->regionlevelmax,1); - /* Read vector attributes*/ - int size; - buf.Read(&size,1); - int* psid2index=xNew(size); - buf.Read(psid2index,size); - this->sid2index.clear(); - this->sid2index.assign(psid2index,psid2index+size); - /* Read geometric mesh (father)*/ - TPZSaveable *sv1 = TPZSaveable::Restore(buf,0); - this->fathermesh = dynamic_cast(sv1); - /* Read geometric mesh (current)*/ - TPZSaveable *sv2 = TPZSaveable::Restore(buf,0); - this->currentmesh = dynamic_cast(sv2); - /* Cleanup*/ - xDelete(psid2index); - } - catch(const std::exception& e) - { - _error_("AdaptiveMeshRefinement::Read: Exception catched!\n"); - } - -} -/*}}}*/ -template class TPZRestoreClass;/*{{{*/ -/*}}}*/ -void AdaptiveMeshRefinement::Write(TPZStream &buf,int withclassid){/*{{{*/ - - try - { - /* Write context (this class) class ID*/ - TPZSaveable::Write(buf,withclassid); - /* Write this class id*/ - int classid = this->ClassId(); - buf.Write(&classid,1); - /* Write simple attributes */ - buf.Write(&this->levelmax,1); - buf.Write(&this->elementswidth,1); - buf.Write(&this->regionlevel1,1); - buf.Write(&this->regionlevelmax,1); - /* Write vector attributes*/ - int size=this->sid2index.size(); - int* psid2index=&this->sid2index[0]; - buf.Write(&size,1);//vector size - buf.Write(psid2index,this->sid2index.size()); - /* Write the geometric mesh*/ - this->fathermesh->Write(buf,this->ClassId()); - this->currentmesh->Write(buf,this->ClassId()); - } - catch(const std::exception& e) - { - _error_("AdaptiveMeshRefinement::Write: Exception catched!\n"); - } + this->index2sid.clear(); + this->specialelementsindex.clear(); } /*}}}*/ @@ -176,27 +97,12 @@ if(numberofelements!=this->sid2index.size()) _error_("Impossible to execute refinement: sid2index.size is not equal to numberofelements!\n"); - /*Execute the refinement.*/ - this->RefinementProcess(amr_verbose,partiallyfloatedelements,masklevelset,deviatorictensorerror,thicknesserror); - - /*Get new geometric mesh in ISSM data structure*/ - this->GetMesh(newnumberofvertices,newnumberofelements,x,y,elementslist); - - /*Verify the new geometry*/ - this->CheckMesh(newnumberofvertices,newnumberofelements,this->elementswidth,x,y,elementslist); - -} -/*}}}*/ -void AdaptiveMeshRefinement::RefinementProcess(bool &amr_verbose,double* partiallyfloatedelements,double* masklevelset, - double* deviatorictensorerror,double* thicknesserror){/*{{{*/ - - if(amr_verbose) _printf_("\n\trefinement process started (level max = " << this->levelmax << ")\n"); - - /*Intermediaries*/ - TPZGeoMesh* nohangingnodesmesh=NULL; + /*Combine the fields*/ double mean_mask = 0; double mean_tauerror = 0; double mean_Herror = 0; - double group_error = 0; - + int* fmask = xNew(numberofelements); + int* ftauerror = xNew(numberofelements); + int* fHerror = xNew(numberofelements); + int* phi = xNew(numberofelements); /*Calculate mean values*/ for(int i=0;isid2index.size();i++){ @@ -208,32 +114,183 @@ mean_tauerror /= this->sid2index.size(); mean_Herror /= this->sid2index.size(); - + /*Filter to thickness*/ + for(int i=0;isid2index.size();i++){ + fmask[i]=0; + if(thicknesserror[i]>mean_Herror) fmask[i]=1; + } + /*Filter to tau*/ + for(int i=0;isid2index.size();i++){ + ftauerror[i]=0; + if(deviatorictensorerror[i]>mean_tauerror) ftauerror[i]=1; + } + /*Sum*/ + for(int i=0;isid2index.size();i++){ + phi[i]=ftauerror[i]+fmask[i]; + } + + /*Execute the refinement.*/ + this->RefinementProcess(amr_verbose,partiallyfloatedelements,masklevelset,deviatorictensorerror,thicknesserror); + + /*Get new geometric mesh in ISSM data structure*/ + this->GetMesh(newnumberofvertices,newnumberofelements,x,y,elementslist); + + std::ofstream file("/home/santos/mesh.vtk"); + TPZVTKGeoMesh::PrintGMeshVTK(this->currentmesh,file); + + /*Verify the new geometry*/ + this->CheckMesh(newnumberofvertices,newnumberofelements,this->elementswidth,x,y,elementslist); + +} +/*}}}*/ +void AdaptiveMeshRefinement::RefinementProcess(bool &amr_verbose,double* partiallyfloatedelements,double* masklevelset,double* deviatorictensorerror,double* thicknesserror){/*{{{*/ + + if(amr_verbose) _printf_("\n\trefinement process started (level max = " << this->levelmax << ")\n"); + + /*Intermediaries*/ + TPZGeoMesh* nohangingnodesmesh=NULL; + + //itapopo + this->step++; + + 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(amr_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->currentmesh->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->currentmesh->DeleteElement(siblings[j],siblings[j]->Index()); + } + } + + /*Now, verify if the father should be refined with uniform pattern (smoother)*/ + this->smooth_frequency=10000000; + if(siblings.size()==3 || this->step%this->smooth_frequency==0){//it keeps the mesh with uniform elements + /*Father has uniform subelements now*/ + TPZVec sons; + father->Divide(sons); + this->smooth_frequency=this->step; + }else{ + specialfatherindex.push_back(father->Index()); + } + if(this->specialelementsindex[i]!=-1) _error_("special element "<currentmesh->BuildConnectivity(); + if(amr_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->currentmesh->Element(this->sid2index[i]); - if(geoel->HasSubElement()) _error_("Impossible to refine: geoel (index) "<sid2index[i]<<" has subelements!\n"); - if(geoel->MaterialId()!=this->GetElemMaterialID()) _error_("Impossible to refine: geoel->MaterialId is not GetElemMaterialID!\n"); - /*Refine*/ - if(thicknesserror[i]>mean_Herror){ + 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()levelmax) geoel->Divide(sons); + if(geoel->Level()levelmax && count==0) geoel->Divide(sons); } - else if(geoel->Level()>0){ /*try to unrefine*/ - TPZVec sons; - geoel->Father()->GetHigherSubElements(sons); - group_error=0; - for(int j=0;jIndex(); + 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->currentmesh->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->currentmesh->DeleteElement(siblings[j],siblings[j]->Index()); + }//for j + }//if + }/*Unrefine process*/ + }//for i this->currentmesh->BuildConnectivity(); if(amr_verbose) _printf_("\t\trefine to avoid hanging nodes...\n"); this->RefineMeshToAvoidHangingNodes(this->currentmesh); - this->currentmesh->BuildConnectivity(); //nohangingnodesmesh = this->CreateRefPatternMesh(newmesh); itapopo tentar otimizar - + if(amr_verbose) _printf_("\trefinement process done!\n"); } @@ -257,5 +314,6 @@ void AdaptiveMeshRefinement::RefineMeshToAvoidHangingNodes(TPZGeoMesh *gmesh){/*{{{*/ - /*Refine elements to avoid hanging nodes: non-uniform refinement*/ + /*Now, insert special elements to avoid hanging nodes*/ + this->specialelementsindex.clear(); const int NElem = gmesh->NElements(); for(int i=0;i refp=TPZRefPatternTools::PerfectMatchRefPattern(geoel); if(refp){ + /*Non-uniform refinement*/ TPZVec Sons; geoel->SetRefPattern(refp); geoel->Divide(Sons); - } - } - gmesh->BuildConnectivity(); - + /*Keep the index of the special elements*/ + for(int j=0;jspecialelementsindex.push_back(Sons[j]->Index()); + } + } + gmesh->BuildConnectivity(); } /*}}}*/ @@ -281,6 +341,6 @@ /* IMPORTANT! pelements are in Matlab indexing NEOPZ works only in C indexing. - This method cleans up and updated the this->sid2index - and fills in it with the new mesh. + This method cleans up and updated the this->sid2index + and this->index2sid and fills in it with the new mesh. Avoid to call this method before Refinement Process.*/ @@ -295,4 +355,5 @@ TPZGeoEl* geoel = NULL; long* vertex_index2sid = xNew(gmesh->NNodes()); + this->index2sid.clear(); this->index2sid.resize(gmesh->NElements()); this->sid2index.clear(); @@ -307,4 +368,7 @@ /*Fill in the vertex_index2sid vector with non usual index value*/ for(int i=0;iNNodes();i++) vertex_index2sid[i]=-1; + + /*Fill in the this->index2sid vector with non usual index value*/ + for(int i=0;iNElements();i++) this->index2sid[i]=-1; /*Get elements without sons and fill in the vertex_index2sid with used vertices (indexes) */ @@ -316,4 +380,5 @@ if(geoel->MaterialId() != this->GetElemMaterialID()) continue; this->sid2index.push_back(i);//keep the element index + this->index2sid[i]=this->sid2index.size()-1;//keep the element sid for(int j=0;jelementswidth;j++){ nodeindex=geoel->NodeIndex(j); @@ -347,6 +412,27 @@ newelements[i*this->elementswidth+j]=(int)sid+1;//C to Matlab indexing } - } - + /*Verify the Jacobian determinant. If detJ<0, swap the 2 first postions: + a -> b + b -> a */ + double detJ,xa,xb,xc,ya,yb,yc; + int a,b,c; + + a=newelements[i*this->elementswidth+0]-1; + b=newelements[i*this->elementswidth+1]-1; + c=newelements[i*this->elementswidth+2]-1; + + xa=newmeshX[a]; ya=newmeshY[a]; + xb=newmeshX[b]; yb=newmeshY[b]; + xc=newmeshX[c]; yc=newmeshY[c]; + + detJ=(xb-xa)*(yc-ya)-(xc-xa)*(yb-ya); + + /*verify and swap, if necessary*/ + if(detJ<0) { + newelements[i*this->elementswidth+0]=b+1;//a->b + newelements[i*this->elementswidth+1]=a+1;//b->a + } + } + /*Setting outputs*/ nvertices = nconformvertices; @@ -358,5 +444,4 @@ /*Cleanup*/ xDelete(vertex_index2sid); - } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h =================================================================== --- /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h (revision 21918) +++ /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h (revision 21919) @@ -12,5 +12,5 @@ /*REAL and STATE definitions, NeoPZ variables itapopo should be read by NeoPZ's config.h*/ #ifndef REFPATTERNDIR - # define REFPATTERNDIR "/home/santos/trunk-jpl/externalpackages/neopz/install/include/refpatterns" + #define REFPATTERNDIR "/home/santos/trunk-jpl/externalpackages/neopz/install/include/refpatterns" #endif @@ -24,5 +24,4 @@ #include -#include #include #include @@ -36,8 +35,13 @@ #include #include +#include #include +#include + +#include "../shared/shared.h" + /*}}}*/ -class AdaptiveMeshRefinement : public TPZSaveable { +class AdaptiveMeshRefinement{ public: @@ -49,9 +53,4 @@ AdaptiveMeshRefinement & operator= (const AdaptiveMeshRefinement &cp); virtual ~AdaptiveMeshRefinement(); - - /*Savable methods*/ - virtual int ClassId() const; - virtual void Read(TPZStream &buf,void *context); - virtual void Write(TPZStream &buf,int withclassid); /*General methods*/ @@ -70,9 +69,13 @@ private: /*Private attributes*/ - int elementswidth; // Geometric nodes for element: 3 == Tria, 4 == Tetra, 6 == Penta + int elementswidth; // Geometric nodes for element: 3 == Tria, 4 == Tetra, 6 == Penta int levelmax; // Max level of refinement double regionlevel1; // Region which will be refined with level 1 double regionlevelmax; // Region which will be refined with level max + int step; //itapopo testando + int smooth_frequency; //itapopo testando 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 *currentmesh; // Current Mesh is the refined mesh Index: /issm/trunk-jpl/src/c/classes/AmrBamg.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/AmrBamg.cpp (revision 21918) +++ /issm/trunk-jpl/src/c/classes/AmrBamg.cpp (revision 21919) @@ -17,10 +17,23 @@ /*Constructor, copy, clean up and destructor*/ -AmrBamg::AmrBamg(IssmDouble hmin, IssmDouble hmax,int fieldenum_in,IssmDouble err){/*{{{*/ +AmrBamg::AmrBamg(IssmDouble hmin,IssmDouble hmax,int fieldenum_in,IssmDouble err_in,int keepmetric_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){/*{{{*/ - this->fieldenum = fieldenum_in; - this->geometry = NULL; - this->fathermesh = NULL; - this->previousmesh = NULL; + 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; + this->geometry = NULL; + this->fathermesh = NULL; + this->previousmesh = NULL; /*Only initialize options for now (same as bamg.m)*/ @@ -49,5 +62,5 @@ this->options->errSize[0]=1; this->options->errSize[1]=1; - this->options->err[0] = err; + this->options->err[0] = err_in; } /*}}}*/ @@ -83,5 +96,5 @@ delete Th; }/*}}}*/ -void AmrBamg::ExecuteRefinementBamg(IssmDouble* field,int* pnewnumberofvertices,int *pnewnumberofelements,IssmDouble** px,IssmDouble** py,IssmDouble** pz,int** pelementslist){/*{{{*/ +void AmrBamg::ExecuteRefinementBamg(IssmDouble* field,IssmDouble* hmaxVertices,int* pnewnumberofvertices,int *pnewnumberofelements,IssmDouble** px,IssmDouble** py,IssmDouble** pz,int** pelementslist){/*{{{*/ /*Intermediaries*/ @@ -93,25 +106,40 @@ _assert_(this->options); _assert_(this->fathermesh); - _assert_(field); + _assert_(field || hmaxVertices);//at least one is necessary /*Prepare field for metric*/ - this->options->field = field; + this->options->field = field; + this->options->hmaxVertices = hmaxVertices; /*remesh*/ if(this->previousmesh){ - this->options->fieldSize[0] = this->previousmesh->VerticesSize[0]; - this->options->fieldSize[1] = 1; + this->options->fieldSize[0] = this->previousmesh->VerticesSize[0]; + this->options->fieldSize[1] = 1; + this->options->hmaxVerticesSize[0] = this->previousmesh->VerticesSize[0]; + this->options->hmaxVerticesSize[1] = 1; Bamgx(meshout,geomout,this->previousmesh,this->geometry,this->options); } else{ - this->options->fieldSize[0] = this->fathermesh->VerticesSize[0]; - this->options->fieldSize[1] = 1; + this->options->fieldSize[0] = this->fathermesh->VerticesSize[0]; + this->options->fieldSize[1] = 1; + this->options->hmaxVerticesSize[0] = this->fathermesh->VerticesSize[0]; + this->options->hmaxVerticesSize[1] = 1; Bamgx(meshout,geomout,this->fathermesh,this->geometry,this->options); } - /*remove field for memory management (FemModel is taking care of deleting it)*/ + /*remove field and hmaxVertices for memory management (FemModel is taking care of deleting it)*/ this->options->field = NULL; this->options->fieldSize[0] = 0; this->options->fieldSize[1] = 0; + this->options->hmaxVertices = NULL; + this->options->hmaxVerticesSize[0] = 0; + this->options->hmaxVerticesSize[1] = 0; + + /*verify if the metric will be reseted or not*/ + if(!this->keepmetric){ + if(this->options->metric) xDelete(this->options->metric); + this->options->metricSize[0] = 0; + this->options->metricSize[1] = 0; + } /*Change previous mesh*/ Index: /issm/trunk-jpl/src/c/classes/AmrBamg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/AmrBamg.h (revision 21918) +++ /issm/trunk-jpl/src/c/classes/AmrBamg.h (revision 21919) @@ -14,12 +14,25 @@ public: int fieldenum; + int keepmetric; + IssmDouble groundingline_resolution; + IssmDouble groundingline_distance; + IssmDouble icefront_resolution; + IssmDouble icefront_distance; + IssmDouble thicknesserror_resolution; + IssmDouble thicknesserror_threshold; + IssmDouble deviatoricerror_resolution; + IssmDouble deviatoricerror_threshold; /* Constructor, destructor etc*/ - AmrBamg(IssmDouble hmin, IssmDouble hmax,int fieldenum_in,IssmDouble err_in); + AmrBamg(IssmDouble hmin,IssmDouble hmax,int fieldenum_in,IssmDouble err_in,int keepmetric_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(); /*General methods*/ void Initialize(int* elements,IssmDouble* x,IssmDouble* y,int numberofvertices,int numberofelements); - void ExecuteRefinementBamg(IssmDouble* field,int* pnewnumberofvertices,int *pnewnumberofelements,IssmDouble** px,IssmDouble** py,IssmDouble** pz,int** pelementslist); + void ExecuteRefinementBamg(IssmDouble* field,IssmDouble* hmaxVertices,int* pnewnumberofvertices,int *pnewnumberofelements,IssmDouble** px,IssmDouble** py,IssmDouble** pz,int** pelementslist); private: Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 21918) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 21919) @@ -1711,5 +1711,6 @@ name==LevelsetfunctionSlopeYEnum || name==LevelsetfunctionPicardEnum || - //name==CalvingCalvingrateEnum || + name==CalvingCalvingrateEnum || + name==CalvingMeltingrateEnum || name==GradientEnum || name==OldGradientEnum || Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 21918) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 21919) @@ -2303,7 +2303,7 @@ #endif - #if defined(_HAVE_BAMG_) && !defined(_HAVE_ADOLC_) + #if defined(_HAVE_BAMG_) && !defined(_HAVE_ADOLC_) case AmrBamgEnum: this->ReMeshBamg(&newnumberofvertices,&newnumberofelements,&newx,&newy,&newz,&newelementslist); break; - #endif + #endif default: _error_("not implemented yet"); @@ -2720,9 +2720,7 @@ void FemModel::WriteMeshInResults(void){/*{{{*/ - //itapopo #if defined(_HAVE_NEOPZ_) && !defined(_HAVE_ADOLC_) this->WriteErrorEstimatorsInResults(); #endif - //itapopo int step = -1; @@ -2754,28 +2752,4 @@ this->results->AddResult(new GenericExternalResult(this->results->Size()+1,MeshYEnum, y,numberofvertices,1,step,time)); - - //itapopo - #if defined(_HAVE_NEOPZ_) && !defined(_HAVE_ADOLC_) - if(IssmComm::GetRank()==0){ - TPZFileStream fstr; - std::stringstream ss; - int frictionlaw,levelmax; - this->parameters->FindParam(&frictionlaw,FrictionLawEnum); - this->parameters->FindParam(&levelmax,AmrLevelMaxEnum); - ss << levelmax; - if(frictionlaw==1){ - ss << "_viscous/amr.txt"; - }else if(frictionlaw==7){ - ss << "_tsai/amr.txt"; - }else{ - _error_("friction law not supported here."); - } - std::string AMRfile = "/home/santos/Misomip2/L" + ss.str(); - fstr.OpenWrite(AMRfile.c_str()); - int withclassid = 1; - this->amr->Write(fstr,withclassid); - } - #endif - //itapopo /*Cleanup*/ @@ -3614,4 +3588,45 @@ } /*}}}*/ +void FemModel::GetElementCenterCoordinates(IssmDouble** pxc,IssmDouble** pyc){/*{{{*/ + + /*Intermediaries*/ + int elementswidth = this->GetElementsWidth(); + int numberofelements = this->elements->NumberOfElements(); + int* elem_vertices = xNew(elementswidth); + 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) ; + + /*Insert the element center coordinates*/ + for(int i=0;ielements->Size();i++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); + 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); + } + + /*Assemble*/ + vxc->Assemble(); + vyc->Assemble(); + + /*Serialize and set output*/ + (*pxc)=vxc->ToMPISerial(); + (*pyc)=vyc->ToMPISerial(); + + /*Cleanup*/ + xDelete(x); + xDelete(y); + xDelete(z); + xDelete(elem_vertices); + delete vxc; + delete vyc; +} +/*}}}*/ #endif @@ -4336,18 +4351,39 @@ int my_rank = IssmComm::GetRank(); + /*Intermediaries*/ + int numberofvertices = this->vertices->NumberOfVertices(); + IssmDouble* vector_serial = NULL; + IssmDouble* hmaxvertices_serial = NULL; + Vector *vector = NULL; + /*Get vector to create metric*/ - int numberofvertices = this->vertices->NumberOfVertices(); - Vector *vector = NULL; - GetVectorFromInputsx(&vector,this,this->amrbamg->fieldenum,VertexSIdEnum); - vector->Assemble(); - IssmDouble* vector_serial = vector->ToMPISerial(); + if(this->amrbamg->fieldenum!=NoneEnum){ + GetVectorFromInputsx(&vector,this,this->amrbamg->fieldenum,VertexSIdEnum); + vector->Assemble(); + vector_serial = vector->ToMPISerial(); + } + + /*Get hmaxVertices to create metric*/ + if(this->amrbamg->groundingline_distance>0||this->amrbamg->icefront_distance>0|| + this->amrbamg->thicknesserror_threshold>0||this->amrbamg->deviatoricerror_threshold>0){ + /*Initialize hmaxvertices with NAN*/ + hmaxvertices_serial=xNew(numberofvertices); + for(int i=0;iamrbamg->groundingline_distance>0) this->GethmaxVerticesFromZeroLevelSetDistance(hmaxvertices_serial,MaskGroundediceLevelsetEnum); + if(this->amrbamg->icefront_distance>0) this->GethmaxVerticesFromZeroLevelSetDistance(hmaxvertices_serial,MaskIceLevelsetEnum); + if(this->amrbamg->thicknesserror_threshold>0) this->GethmaxVerticesFromEstimators(hmaxvertices_serial,ThicknessErrorEstimatorEnum); + 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); + if(newnumberofvertices<=0 || newnumberofelements<=0) _error_("Error in the refinement process."); + } + + /*Cleanup*/ + xDelete(vector_serial); + xDelete(hmaxvertices_serial); delete vector; - - if(my_rank==0){ - this->amrbamg->ExecuteRefinementBamg(vector_serial,&newnumberofvertices,&newnumberofelements,&newx,&newy,&newz,&newelementslist); - if(newnumberofvertices<=0 || newnumberofelements<=0) _error_("Error in the refinement process."); - } - - xDelete(vector_serial); /*Send new mesh to others CPU*/ @@ -4379,12 +4415,12 @@ int numberofvertices = this->vertices->NumberOfVertices(); int numberofelements = this->elements->NumberOfElements(); - int numberofsegments = 0; //used on matlab IssmDouble* x = NULL; IssmDouble* y = NULL; IssmDouble* z = NULL; int* elements = NULL; - int elementswidth = this->GetElementsWidth(); //just tria elements in this version. Itapopo: IssmDouble hmin,hmax,err; - int fieldenum; + IssmDouble groundingline_resolution,groundingline_distance,icefront_resolution,icefront_distance; + IssmDouble thicknesserror_resolution,thicknesserror_threshold,deviatoricerror_resolution,deviatoricerror_threshold; + int fieldenum,keepmetric; /*Get rank*/ @@ -4402,7 +4438,20 @@ this->parameters->FindParam(&fieldenum,AmrFieldEnum); this->parameters->FindParam(&err,AmrErrEnum); + this->parameters->FindParam(&keepmetric,AmrKeepMetricEnum); + 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); + this->amrbamg = new AmrBamg(hmin,hmax,fieldenum,err,keepmetric, + groundingline_resolution,groundingline_distance, + icefront_resolution,icefront_distance, + thicknesserror_resolution,thicknesserror_threshold, + deviatoricerror_resolution,deviatoricerror_threshold); /*Re-create original mesh and put it in bamg structure (only cpu 0)*/ @@ -4416,4 +4465,202 @@ xDelete(z); xDelete(elements); +} +/*}}}*/ +void FemModel::GethmaxVerticesFromZeroLevelSetDistance(IssmDouble* hmaxvertices,int levelset_type){/*{{{*/ + + if(!hmaxvertices) _error_("hmaxvertices is NULL!\n"); + + /*Intermediaries*/ + int numberofvertices = this->vertices->NumberOfVertices(); + IssmDouble* verticedistance = NULL; + IssmDouble threshold,resolution; + + switch(levelset_type){ + case MaskGroundediceLevelsetEnum: + threshold = this->amrbamg->groundingline_distance; + resolution = this->amrbamg->groundingline_resolution; + break; + case MaskIceLevelsetEnum: + threshold = this->amrbamg->icefront_distance; + resolution = this->amrbamg->icefront_resolution; + break; + default: _error_("not implemented yet"); + } + + /*Get vertice distance to zero level set points*/ + this->GetVerticeDistanceToZeroLevelSet(&verticedistance,levelset_type); + if(!verticedistance) _error_("verticedistance is NULL!\n"); + + /*Fill hmaxVertices*/ + for(int i=0;i(hmaxvertices[i])) hmaxvertices[i]=resolution; + else hmaxvertices[i]=min(resolution,hmaxvertices[i]); + } + } + + /*Cleanup*/ + xDelete(verticedistance); +} +/*}}}*/ +void FemModel::GethmaxVerticesFromEstimators(IssmDouble* hmaxvertices,int errorestimator_type){/*{{{*/ + + if(!hmaxvertices) _error_("hmaxvertices is NULL!\n"); + + /*Intermediaries*/ + int numberofelements = this->elements->NumberOfElements(); + int numberofvertices = this->vertices->NumberOfVertices(); + IssmDouble* error_elements = NULL; + IssmDouble* error_vertices = NULL; + IssmDouble *x = NULL; + IssmDouble *y = NULL; + IssmDouble *z = NULL; + int *index = NULL; + IssmDouble maxerror,threshold,resolution; + + switch(errorestimator_type){ + case ThicknessErrorEstimatorEnum: + threshold = this->amrbamg->thicknesserror_threshold; + resolution = this->amrbamg->thicknesserror_resolution; + this->ThicknessZZErrorEstimator(&error_elements); + break; + case DeviatoricStressErrorEstimatorEnum: + threshold = this->amrbamg->deviatoricerror_threshold; + resolution = this->amrbamg->deviatoricerror_resolution; + this->ZZErrorEstimator(&error_elements); + break; + default: _error_("not implemented yet"); + } + + if(!error_elements) _error_("error_elements is NULL!\n"); + + /*Get mesh*/ + this->GetMesh(this->vertices,this->elements,&x,&y,&z,&index); + + /*Sum the estimators over the vertices*/ + error_vertices=xNewZeroInit(numberofvertices); + for(int i=0;iGetElementsWidth();j++){ + int vid=index[i*this->GetElementsWidth()+j]-1;//Matlab to C indexing + error_vertices[vid]+=error_elements[i]; + } + } + + /*Find the max of the estimators (use error_elements)*/ + maxerror=error_elements[0]; + for(int i=0;ithreshold*maxerror){ + if(xIsNan(hmaxvertices[i])) hmaxvertices[i]=resolution; + else hmaxvertices[i]=min(resolution,hmaxvertices[i]); + } + } + + /*Cleanup*/ + xDelete(x); + xDelete(y); + xDelete(z); + xDelete(index); + xDelete(error_elements); + xDelete(error_vertices); +} +/*}}}*/ +void FemModel::GetVerticeDistanceToZeroLevelSet(IssmDouble** pverticedistance,int levelset_type){/*{{{*/ + + /*Here, "zero level set" means grounding line or ice front, depending on the level set type*/ + /*pverticedistance is the minimal vertice distance to the grounding line or ice front*/ + if(levelset_type!=MaskGroundediceLevelsetEnum && levelset_type!=MaskIceLevelsetEnum){ + _error_("level set type not implemented yet!"); + } + + /*Output*/ + IssmDouble* verticedistance; + + /*Intermediaries*/ + int elementswidth = this->GetElementsWidth(); + int numberofvertices = this->vertices->NumberOfVertices(); + int numberofelements = this->elements->NumberOfElements(); + int* elem_vertices = xNew(elementswidth); + IssmDouble* levelset = xNew(elementswidth); + IssmDouble* xp = NULL; + IssmDouble* yp = NULL; + IssmDouble* x = NULL; + IssmDouble* y = NULL; + IssmDouble* z = NULL; + Vector* vx_zerolevelset = new Vector(numberofelements); + Vector* vy_zerolevelset = new Vector(numberofelements); + IssmDouble* x_zerolevelset = NULL; + IssmDouble* y_zerolevelset = NULL; + int count,sid,npoints; + IssmDouble xcenter,ycenter,distance; + + /*Get vertices coordinates*/ + VertexCoordinatesx(&x,&y,&z,this->vertices,false) ; + + /*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(); + xcenter = NAN; + ycenter = 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); + } + /*Assemble and serialize*/ + vx_zerolevelset->Assemble(); + vy_zerolevelset->Assemble(); + x_zerolevelset=vx_zerolevelset->ToMPISerial(); + y_zerolevelset=vy_zerolevelset->ToMPISerial(); + + /*keep just the element center coordinates with zero level set (compact the structure to save time in the next step)*/ + count = 0; + xp = xNewZeroInit(numberofelements); + yp = xNewZeroInit(numberofelements); + for(int i=0;i(x_zerolevelset[i])){ + xp[count]=x_zerolevelset[i]; + yp[count]=y_zerolevelset[i]; + count++; + } + } + npoints=count; + + /*Find the minimal vertice distance to the zero levelset (grounding line or ice front)*/ + verticedistance=xNew(numberofvertices); + for(int i=0;i(elem_vertices); + xDelete(levelset); + xDelete(x_zerolevelset); + xDelete(y_zerolevelset); + xDelete(xp); + xDelete(yp); + xDelete(x); + xDelete(y); + xDelete(z); + delete vx_zerolevelset; + delete vy_zerolevelset; } /*}}}*/ @@ -4480,5 +4727,4 @@ xDelete(deviatorictensorerror); xDelete(thicknesserror); - } /*}}}*/ @@ -4512,35 +4758,9 @@ /*Create initial mesh (coarse mesh) in neopz data structure*/ /*Just CPU #0 should keep AMR object*/ - this->SetRefPatterns(); + /*Initialize refinement pattern*/ + this->SetRefPatterns(); if(my_rank==0){ - bool ismisomip = false; - if(ismisomip){//itapopo - TPZFileStream fstr; - std::stringstream ss; - int frictionlaw; - - this->parameters->FindParam(&frictionlaw,FrictionLawEnum); - - ss << levelmax; - if(frictionlaw==1){ - ss << "_viscous/amr.txt"; - }else if(frictionlaw==7){ - ss << "_tsai/amr.txt"; - }else{ - _error_("friction law not supported here."); - } - - std::string AMRfile = "/home/santos/Misomip2/L" + ss.str(); - fstr.OpenRead(AMRfile.c_str()); - - TPZSaveable *sv = TPZSaveable::Restore(fstr,0); - this->amr = dynamic_cast(sv); - } - else{ - this->amr = new AdaptiveMeshRefinement(); - //this->amr->SetLevelMax(levelmax); //Set max level of refinement - //this->amr->SetRegions(regionlevel1,regionlevelmax); - this->amr->CreateInitialMesh(numberofvertices,numberofelements,elementswidth,x,y,elements); - } + this->amr = new AdaptiveMeshRefinement(); + this->amr->CreateInitialMesh(numberofvertices,numberofelements,elementswidth,x,y,elements); this->amr->SetLevelMax(levelmax); //Set max level of refinement this->amr->SetRegions(regionlevel1,regionlevelmax); @@ -4559,13 +4779,15 @@ gRefDBase.InitializeUniformRefPattern(ETriangle); - /*Insert specifics patterns to ISSM core*/ + /*Insert specifics patterns to ISSM core*/ std::string filepath = REFPATTERNDIR; std::string filename1 = filepath + "/2D_Triang_Rib_3.rpt"; std::string filename2 = filepath + "/2D_Triang_Rib_4.rpt"; - std::string filename3 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_4.rpt"; - std::string filename4 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_4_permuted.rpt"; - std::string filename5 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_5.rpt"; - std::string filename6 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_5_permuted.rpt"; - std::string filename7 = filepath + "/2D_Triang_Rib_5.rpt"; + std::string filename3 = filepath + "/2D_Triang_Rib_5.rpt"; + std::string filename4 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_4.rpt"; + std::string filename5 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_4_permuted.rpt"; + std::string filename6 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_5.rpt"; + std::string filename7 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_5_permuted.rpt"; + std::string filename8 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_4_5.rpt"; + std::string filename9 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_4_5_permuted.rpt"; TPZAutoPointer refpat1 = new TPZRefPattern(filename1); @@ -4576,4 +4798,6 @@ TPZAutoPointer refpat6 = new TPZRefPattern(filename6); TPZAutoPointer refpat7 = new TPZRefPattern(filename7); + TPZAutoPointer refpat8 = new TPZRefPattern(filename8); + TPZAutoPointer refpat9 = new TPZRefPattern(filename9); if(!gRefDBase.FindRefPattern(refpat1)) gRefDBase.InsertRefPattern(refpat1); @@ -4584,4 +4808,6 @@ if(!gRefDBase.FindRefPattern(refpat6)) gRefDBase.InsertRefPattern(refpat6); if(!gRefDBase.FindRefPattern(refpat7)) gRefDBase.InsertRefPattern(refpat7); + if(!gRefDBase.FindRefPattern(refpat8)) gRefDBase.InsertRefPattern(refpat8); + if(!gRefDBase.FindRefPattern(refpat9)) gRefDBase.InsertRefPattern(refpat9); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 21918) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 21919) @@ -177,4 +177,5 @@ void ThicknessZZErrorEstimator(IssmDouble** pelementerror); void MeanGroundedIceLevelSet(IssmDouble** pmasklevelset); + void GetElementCenterCoordinates(IssmDouble** pxc,IssmDouble** pyc); #endif @@ -182,8 +183,10 @@ void ReMeshBamg(int* pnewnumberofvertices,int* pnewnumberofelements,IssmDouble** pnewx,IssmDouble** pnewy,IssmDouble** pnewz,int** pnewelementslist); void InitializeAdaptiveRefinementBamg(void); + void GethmaxVerticesFromZeroLevelSetDistance(IssmDouble* hmaxvertices,int levelset_type); + void GethmaxVerticesFromEstimators(IssmDouble* hmaxvertices,int errorestimator_type); + void GetVerticeDistanceToZeroLevelSet(IssmDouble** pverticedistance,int leveset_type); #endif #if defined(_HAVE_NEOPZ_) && !defined(_HAVE_ADOLC_) - /*Adaptive mesh refinement methods*/ void ReMeshNeopz(int* pnewnumberofvertices,int* pnewnumberofelements,IssmDouble** pnewx,IssmDouble** pnewy,IssmDouble** pnewz,int** pnewelementslist); void InitializeAdaptiveRefinementNeopz(void); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 21918) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 21919) @@ -114,4 +114,13 @@ parameters->AddObject(iomodel->CopyConstantObject("md.amr.hmax",AmrHmaxEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.amr.err",AmrErrEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.keepmetric",AmrKeepMetricEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.groundingline_resolution",AmrGroundingLineResolutionEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.groundingline_distance",AmrGroundingLineDistanceEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.icefront_resolution",AmrIceFrontResolutionEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.icefront_distance",AmrIceFrontDistanceEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.thicknesserror_resolution",AmrThicknessErrorResolutionEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.thicknesserror_threshold",AmrThicknessErrorThresholdEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.deviatoricerror_resolution",AmrDeviatoricErrorResolutionEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.deviatoricerror_threshold",AmrDeviatoricErrorThresholdEnum)); /*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 21918) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21919) @@ -849,4 +849,13 @@ AmrFieldEnum, AmrErrEnum, + AmrKeepMetricEnum, + AmrGroundingLineResolutionEnum, + AmrGroundingLineDistanceEnum, + AmrIceFrontResolutionEnum, + AmrIceFrontDistanceEnum, + AmrThicknessErrorResolutionEnum, + AmrThicknessErrorThresholdEnum, + AmrDeviatoricErrorResolutionEnum, + AmrDeviatoricErrorThresholdEnum, DeviatoricStressErrorEstimatorEnum, ThicknessErrorEstimatorEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21918) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21919) @@ -825,4 +825,13 @@ case AmrFieldEnum : return "AmrField"; case AmrErrEnum : return "AmrErr"; + case AmrKeepMetricEnum : return "AmrKeepMetric"; + case AmrGroundingLineResolutionEnum : return "AmrGroundingLineResolution"; + case AmrGroundingLineDistanceEnum : return "AmrGroundingLineDistance"; + case AmrIceFrontResolutionEnum : return "AmrIceFrontResolution"; + case AmrIceFrontDistanceEnum : return "AmrIceFrontDistance"; + case AmrThicknessErrorResolutionEnum : return "AmrThicknessErrorResolution"; + case AmrThicknessErrorThresholdEnum : return "AmrThicknessErrorThreshold"; + case AmrDeviatoricErrorResolutionEnum : return "AmrDeviatoricErrorResolution"; + case AmrDeviatoricErrorThresholdEnum : return "AmrDeviatoricErrorThreshold"; case DeviatoricStressErrorEstimatorEnum : return "DeviatoricStressErrorEstimator"; case ThicknessErrorEstimatorEnum : return "ThicknessErrorEstimator"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21918) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21919) @@ -843,4 +843,13 @@ else if (strcmp(name,"AmrField")==0) return AmrFieldEnum; else if (strcmp(name,"AmrErr")==0) return AmrErrEnum; + else if (strcmp(name,"AmrKeepMetric")==0) return AmrKeepMetricEnum; + else if (strcmp(name,"AmrGroundingLineResolution")==0) return AmrGroundingLineResolutionEnum; + else if (strcmp(name,"AmrGroundingLineDistance")==0) return AmrGroundingLineDistanceEnum; + else if (strcmp(name,"AmrIceFrontResolution")==0) return AmrIceFrontResolutionEnum; + else if (strcmp(name,"AmrIceFrontDistance")==0) return AmrIceFrontDistanceEnum; + else if (strcmp(name,"AmrThicknessErrorResolution")==0) return AmrThicknessErrorResolutionEnum; + else if (strcmp(name,"AmrThicknessErrorThreshold")==0) return AmrThicknessErrorThresholdEnum; + else if (strcmp(name,"AmrDeviatoricErrorResolution")==0) return AmrDeviatoricErrorResolutionEnum; + else if (strcmp(name,"AmrDeviatoricErrorThreshold")==0) return AmrDeviatoricErrorThresholdEnum; else if (strcmp(name,"DeviatoricStressErrorEstimator")==0) return DeviatoricStressErrorEstimatorEnum; else if (strcmp(name,"ThicknessErrorEstimator")==0) return ThicknessErrorEstimatorEnum; @@ -866,5 +875,8 @@ else if (strcmp(name,"ExternalResult")==0) return ExternalResultEnum; else if (strcmp(name,"FileParam")==0) return FileParamEnum; - else if (strcmp(name,"Input")==0) return InputEnum; + else stage=8; + } + if(stage==8){ + if (strcmp(name,"Input")==0) return InputEnum; else if (strcmp(name,"IntInput")==0) return IntInputEnum; else if (strcmp(name,"IntParam")==0) return IntParamEnum; @@ -875,8 +887,5 @@ else if (strcmp(name,"Matenhancedice")==0) return MatenhancediceEnum; else if (strcmp(name,"Matestar")==0) return MatestarEnum; - else stage=8; - } - if(stage==8){ - if (strcmp(name,"Matpar")==0) return MatparEnum; + else if (strcmp(name,"Matpar")==0) return MatparEnum; else if (strcmp(name,"Node")==0) return NodeEnum; else if (strcmp(name,"Numericalflux")==0) return NumericalfluxEnum; @@ -989,5 +998,8 @@ else if (strcmp(name,"MaxVel")==0) return MaxVelEnum; else if (strcmp(name,"MinVx")==0) return MinVxEnum; - else if (strcmp(name,"MaxVx")==0) return MaxVxEnum; + else stage=9; + } + if(stage==9){ + if (strcmp(name,"MaxVx")==0) return MaxVxEnum; else if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum; else if (strcmp(name,"MinVy")==0) return MinVyEnum; @@ -998,8 +1010,5 @@ else if (strcmp(name,"MaxAbsVz")==0) return MaxAbsVzEnum; else if (strcmp(name,"FloatingArea")==0) return FloatingAreaEnum; - else stage=9; - } - if(stage==9){ - if (strcmp(name,"GroundedArea")==0) return GroundedAreaEnum; + else if (strcmp(name,"GroundedArea")==0) return GroundedAreaEnum; else if (strcmp(name,"IceMass")==0) return IceMassEnum; else if (strcmp(name,"IceVolume")==0) return IceVolumeEnum; Index: /issm/trunk-jpl/src/m/classes/amrbamg.m =================================================================== --- /issm/trunk-jpl/src/m/classes/amrbamg.m (revision 21918) +++ /issm/trunk-jpl/src/m/classes/amrbamg.m (revision 21919) @@ -10,4 +10,13 @@ fieldname = ''; err = 0.; + keepmetric = 0; + groundingline_resolution = 0.; + groundingline_distance = 0.; + icefront_resolution = 0.; + icefront_distance = 0.; + thicknesserror_resolution = 0.; + thicknesserror_threshold = 0.; + deviatoricerror_resolution = 0.; + deviatoricerror_threshold = 0.; end methods @@ -30,4 +39,17 @@ self.err=3.; + %keep metric? + self.keepmetric=1; + + %other criterias + self.groundingline_resolution=500.; + self.groundingline_distance=0.; + self.icefront_resolution=500.; + self.icefront_distance=0.; + self.thicknesserror_resolution=500.; + self.thicknesserror_threshold=0.; + self.deviatoricerror_resolution=500.; + self.deviatoricerror_threshold=0.; + end % }}} function md = checkconsistency(self,md,solution,analyses) % {{{ @@ -36,4 +58,13 @@ md = checkfield(md,'fieldname','amr.hmin','numel',[1],'>',0,'<',self.hmax,'NaN',1); %md = checkfield(md,'fieldname','amr.fieldname','string',[1]); + md = checkfield(md,'fieldname','amr.keepmetric','numel',[1],'>=',0,'<=',1,'NaN',1); + md = checkfield(md,'fieldname','amr.groundingline_resolution','numel',[1],'>',0,'<',self.hmax,'NaN',1); + md = checkfield(md,'fieldname','amr.groundingline_distance','numel',[1],'>=',0,'NaN',1,'Inf',1); + md = checkfield(md,'fieldname','amr.icefront_resolution','numel',[1],'>',0,'<',self.hmax,'NaN',1); + md = checkfield(md,'fieldname','amr.icefront_distance','numel',[1],'>=',0,'NaN',1,'Inf',1); + 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.deviatoricerror_resolution','numel',[1],'>',0,'<',self.hmax,'NaN',1); + md = checkfield(md,'fieldname','amr.deviatoricerror_threshold','numel',[1],'>=',0,'<=',1,'NaN',1); end % }}} function disp(self) % {{{ @@ -43,4 +74,13 @@ fielddisplay(self,'hmax',['maximum element length']); fielddisplay(self,'fieldname',['name of input that will be used to compute the metric (should be an input of FemModel)']); + fielddisplay(self,'keepmetric',['indicates whether the metric should be kept every remeshing time']); + fielddisplay(self,'groundingline_resolution',['element length near the grounding line']); + fielddisplay(self,'groundingline_distance',['distance around the grounding line which elements will be refined']); + fielddisplay(self,'icefront_resolution',['element length near the ice front']); + fielddisplay(self,'icefront_distance',['distance around the ice front which elements will be refined']); + fielddisplay(self,'thicknesserror_resolution',['element length when thickness error estimator is used']); + fielddisplay(self,'thicknesserror_threshold',['maximum threshold 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']); end % }}} @@ -51,5 +91,14 @@ WriteData(fid,prefix,'object',self,'class','amr','fieldname','hmax','format','Double'); WriteData(fid,prefix,'object',self,'class','amr','fieldname','fieldname','format','String'); - WriteData(fid,prefix,'object',self,'class','err','fieldname','err','format','Double'); + WriteData(fid,prefix,'object',self,'class','amr','fieldname','err','format','Double'); + WriteData(fid,prefix,'object',self,'class','amr','fieldname','keepmetric','format','Integer'); + WriteData(fid,prefix,'object',self,'class','amr','fieldname','groundingline_resolution','format','Double'); + WriteData(fid,prefix,'object',self,'class','amr','fieldname','groundingline_distance','format','Double'); + WriteData(fid,prefix,'object',self,'class','amr','fieldname','icefront_resolution','format','Double'); + WriteData(fid,prefix,'object',self,'class','amr','fieldname','icefront_distance','format','Double'); + 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','deviatoricerror_resolution','format','Double'); + WriteData(fid,prefix,'object',self,'class','amr','fieldname','deviatoricerror_threshold','format','Double'); end % }}} Index: /issm/trunk-jpl/src/m/mesh/bamg.m =================================================================== --- /issm/trunk-jpl/src/m/mesh/bamg.m (revision 21918) +++ /issm/trunk-jpl/src/m/mesh/bamg.m (revision 21919) @@ -25,7 +25,7 @@ % - maxsubdiv : maximum subdivision of exisiting elements (default is 10) % - metric : matrix (numberofnodes x 3) used as a metric -% - Metrictype : 1 -> absolute error c/(err coeff^2) * Abs(H) (default) -% 2 -> relative error c/(err coeff^2) * Abs(H)/max(s,cutoff*max(s)) -% 3 -> rescaled absolute error c/(err coeff^2) * Abs(H)/(smax-smin) +% - Metrictype : 0 -> absolute error c/(err coeff^2) * Abs(H) (default) +% 1 -> relative error c/(err coeff^2) * Abs(H)/max(s,cutoff*max(s)) +% 2 -> rescaled absolute error c/(err coeff^2) * Abs(H)/(smax-smin) % - nbjacoby : correction used by Hessiantype=1 (default is 1) % - nbsmooth : number of metric smoothing procedure (default is 3)