Index: /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp (revision 22099) +++ /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp (revision 22100) @@ -28,10 +28,14 @@ this->CleanUp(); /*Copy all data*/ - this->fathermesh = new TPZGeoMesh(*cp.fathermesh); - this->currentmesh = new TPZGeoMesh(*cp.currentmesh); - this->levelmax = cp.levelmax; - this->elementswidth = cp.elementswidth; - this->regionlevel1 = cp.regionlevel1; - this->regionlevelmax = cp.regionlevelmax; + this->fathermesh = new TPZGeoMesh(*cp.fathermesh); + this->previousmesh = new TPZGeoMesh(*cp.previousmesh); + this->level_max = cp.level_max; + this->radius_level_max = cp.radius_level_max; + this->gradation = cp.gradation; + this->lag = cp.lag; + this->groundingline_distance = cp.groundingline_distance; + this->icefront_distance = cp.icefront_distance; + this->thicknesserror_threshold = cp.thicknesserror_threshold; + this->deviatoricerror_threshold = cp.deviatoricerror_threshold; this->sid2index.clear(); this->sid2index.resize(cp.sid2index.size()); @@ -56,9 +60,13 @@ /*Verify and delete all data*/ if(this->fathermesh) delete this->fathermesh; - if(this->currentmesh) delete this->currentmesh; - this->levelmax = -1; - this->elementswidth = -1; - this->regionlevel1 = -1; - this->regionlevelmax = -1; + if(this->previousmesh) delete this->previousmesh; + this->level_max = -1; + this->radius_level_max = -1; + this->gradation = -1; + this->lag = -1; + this->groundingline_distance = -1; + this->icefront_distance = -1; + this->thicknesserror_threshold = -1; + this->deviatoricerror_threshold = -1; this->sid2index.clear(); this->index2sid.clear(); @@ -69,12 +77,14 @@ /*Set pointers to NULL*/ - this->fathermesh = NULL; - this->currentmesh = NULL; - this->levelmax = -1; - this->elementswidth = -1; - this->regionlevel1 = -1; - this->regionlevelmax = -1; - this->step = 0; - this->smooth_frequency = 1; + this->fathermesh = NULL; + this->previousmesh = NULL; + this->level_max = -1; + this->radius_level_max = -1; + this->gradation = -1; + this->lag = -1; + this->groundingline_distance = -1; + this->icefront_distance = -1; + this->thicknesserror_threshold = -1; + this->deviatoricerror_threshold = -1; this->sid2index.clear(); this->index2sid.clear(); @@ -84,72 +94,50 @@ /*Mesh refinement methods*/ -void AdaptiveMeshRefinement::Execute(bool &amr_verbose, - int &numberofelements, - double* partiallyfloatedelements, - double *masklevelset, - double* deviatorictensorerror, - double* thicknesserror, - int &newnumberofvertices,int &newnumberofelements,double** x,double** y,int** elementslist){/*{{{*/ - - /*IMPORTANT! newelements are in Matlab indexing*/ +void AdaptiveMeshRefinement::ExecuteRefinement(int numberofpoints,double* xylist,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->currentmesh) _error_("Impossible to execute refinement: fathermesh or currentmesh is NULL!\n"); - if(numberofelements!=this->sid2index.size()) _error_("Impossible to execute refinement: sid2index.size is not equal to numberofelements!\n"); - - /*Combine the fields*/ - double mean_mask = 0; - double mean_tauerror = 0; - double mean_Herror = 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++){ - 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(); - /*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); + if(!this->fathermesh || !this->previousmesh) _error_("Impossible to execute refinement: fathermesh or previousmesh 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); /*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); + this->GetMesh(this->previousmesh,pnewnumberofvertices,pnewnumberofelements,px,py,pelementslist); /*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){/*{{{*/ + this->CheckMesh(pnewnumberofvertices,pnewnumberofelements,px,py,pelementslist); +} +/*}}}*/ +void AdaptiveMeshRefinement::RefinementProcess(bool &verbose,double* gl_elementdistance,double* if_elementdistance,double* deviatoricerror,double* thicknesserror){/*{{{*/ - if(amr_verbose) _printf_("\n\trefinement process started (level max = " << this->levelmax << ")\n"); + if(verbose) _printf_("\n\trefinement process started (level max = " << this->level_max << ")\n"); /*Intermediaries*/ TPZGeoMesh* nohangingnodesmesh=NULL; - - //itapopo - this->step++; double mean_mask = 0; @@ -161,4 +149,5 @@ /*Calculate mean values*/ + /* for(int i=0;isid2index.size();i++){ mean_mask += masklevelset[i]; @@ -169,11 +158,11 @@ mean_tauerror /= this->sid2index.size(); mean_Herror /= this->sid2index.size(); - - if(amr_verbose) _printf_("\t\tdeal with special elements...\n"); + */ + 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->currentmesh->Element(this->specialelementsindex[i]); + TPZGeoEl* geoel=this->previousmesh->Element(this->specialelementsindex[i]); if(!geoel)_error_("special element (sid) "<HasSubElement())_error_("special element (sid) "<NSubElements()<<" subelements!\n"); @@ -213,15 +202,13 @@ /*Ok, the special element can be deleted*/ siblings[j]->ResetSubElements(); - this->currentmesh->DeleteElement(siblings[j],siblings[j]->Index()); + this->previousmesh->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 + if(siblings.size()==3){//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()); @@ -229,12 +216,12 @@ if(this->specialelementsindex[i]!=-1) _error_("special element "<currentmesh->BuildConnectivity(); - - if(amr_verbose) _printf_("\t\tuniform refinement...\n"); + this->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->currentmesh->Element(this->sid2index[i]); + 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"); @@ -252,5 +239,5 @@ } TPZVec sons; - if(geoel->Level()levelmax && count==0) geoel->Divide(sons); + if(geoel->Level()level_max && count==0) geoel->Divide(sons); } else if(geoel->Level()>0) @@ -274,5 +261,5 @@ if(group_Herror>0 && group_Herror<0*mean_Herror){ //itapopo /*Reset subelements in the father*/ - this->currentmesh->Element(geoel->Father()->Index())->ResetSubElements(); + this->previousmesh->Element(geoel->Father()->Index())->ResetSubElements(); /*Delete the elements and set their indexes in the index2sid and sid2index*/ for (int j=0;jindex2sid[index]=-1; if(sid!=-1) this->sid2index[sid]=-1; - this->currentmesh->DeleteElement(siblings[j],siblings[j]->Index()); + this->previousmesh->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->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(amr_verbose) _printf_("\trefinement process done!\n"); -} -/*}}}*/ -void AdaptiveMeshRefinement::RefineMesh(TPZGeoMesh *gmesh,std::vector &elements){/*{{{*/ + 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.;; + TPZVec qsi(2,0.),cp(3,0.); + TPZVec sons; + + /*First, delete the special elements*/ + this->DeleteSpecialElements(verbose,gmesh); + + /*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 "<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)->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*/ @@ -312,30 +403,63 @@ } /*}}}*/ -void AdaptiveMeshRefinement::RefineMeshToAvoidHangingNodes(TPZGeoMesh *gmesh){/*{{{*/ +void AdaptiveMeshRefinement::RefineMeshToAvoidHangingNodes(bool &verbose,TPZGeoMesh* gmesh){/*{{{*/ /*Now, insert special elements to avoid hanging nodes*/ + if(verbose) _printf_("\trefining to avoid hanging nodes (total: "); + + /*Intermediaries*/ + int nelem=-1; + int count= 1; + + while(count>0){ + nelem=gmesh->NElements();//must keep here + count=0; + for(int i=0;iElement(i); + if(!geoel) continue; + if(geoel->HasSubElement()) continue; + if(geoel->MaterialId() != this->GetElemMaterialID()) continue; + /*Get the refinement pattern for this element and refine it*/ + TPZAutoPointer refp=TPZRefPatternTools::PerfectMatchRefPattern(geoel); + if(refp){ + /*Non-uniform refinement*/ + TPZVec sons; + geoel->SetRefPattern(refp); + geoel->Divide(sons); + count++; + /*Keep the index of the special elements*/ + for(int j=0;jspecialelementsindex.push_back(sons[j]->Index()); + } + } + if(verbose){ + if(count==0) _printf_(""<BuildConnectivity(); + } +} +/*}}}*/ +void AdaptiveMeshRefinement::DeleteSpecialElements(bool &verbose,TPZGeoMesh* gmesh){/*{{{*/ + + /*Intermediaries*/ + int count=0; + + if(verbose) _printf_("\tdelete "<specialelementsindex.size()<<" special elements (total: "); + for(int i=0;ispecialelementsindex.size();i++){ + if(this->specialelementsindex[i]==-1) continue; + if(!gmesh->Element(this->specialelementsindex[i])) continue; + if(gmesh->Element(this->specialelementsindex[i])->Father()) gmesh->Element(this->specialelementsindex[i])->Father()->ResetSubElements(); + gmesh->DeleteElement(gmesh->Element(this->specialelementsindex[i]),this->specialelementsindex[i]); + count++; + } + if(verbose) _printf_(""<specialelementsindex.clear(); - const int NElem = gmesh->NElements(); - for(int i=0;iElement(i); - if(!geoel) continue; - if(geoel->HasSubElement()) continue; - if(geoel->MaterialId() != this->GetElemMaterialID()) continue; - /*Get the refinement pattern for this element and refine it*/ - TPZAutoPointer refp=TPZRefPatternTools::PerfectMatchRefPattern(geoel); - if(refp){ - /*Non-uniform refinement*/ - TPZVec Sons; - geoel->SetRefPattern(refp); - geoel->Divide(Sons); - /*Keep the index of the special elements*/ - for(int j=0;jspecialelementsindex.push_back(Sons[j]->Index()); - } - } + /*Adjust connectivities*/ gmesh->BuildConnectivity(); } /*}}}*/ -void AdaptiveMeshRefinement::GetMesh(int &nvertices,int &nelements,double** px,double** py, int** pelements){/*{{{*/ +void AdaptiveMeshRefinement::GetMesh(TPZGeoMesh* gmesh,int* nvertices,int* nelements,double** px,double** py, int** pelements){/*{{{*/ /* IMPORTANT! pelements are in Matlab indexing @@ -346,24 +470,15 @@ /*Intermediaries */ - TPZGeoMesh* gmesh = this->currentmesh;//itapopo confirmar - long sid,nodeindex; - int nconformelements,nconformvertices; + int nconformelements,nconformvertices; + int ntotalvertices = gmesh->NNodes(); int* newelements = NULL; - double* newmeshX = NULL;//xNew(ntotalvertices); - double* newmeshY = NULL;//xNew(ntotalvertices); + double* newmeshX = NULL; + double* newmeshY = NULL; TPZGeoEl* geoel = NULL; - long* vertex_index2sid = xNew(gmesh->NNodes()); + long* vertex_index2sid = xNew(ntotalvertices); this->index2sid.clear(); this->index2sid.resize(gmesh->NElements()); this->sid2index.clear(); - /*Get mesh coords */ - //for(int i=0;i coords(3,0.); - // gmesh->NodeVec()[i].GetCoordinates(coords); - // newmeshX[i] = coords[0]; - // newmeshY[i] = coords[1]; - //} - /*Fill in the vertex_index2sid vector with non usual index value*/ for(int i=0;iNNodes();i++) vertex_index2sid[i]=-1; @@ -379,9 +494,10 @@ if(geoel->HasSubElement()) continue; 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++){ + this->sid2index.push_back(geoel->Index());//keep the element index + this->index2sid[geoel->Index()]=this->sid2index.size()-1;//keep the element sid + for(int j=0;jGetNumberOfNodes();j++){ nodeindex=geoel->NodeIndex(j); - if(vertex_index2sid[nodeindex]==-1){ + if(nodeindex<0) _error_("nodeindex is <0\n"); + if(vertex_index2sid[nodeindex]==-1){ vertex_index2sid[nodeindex]=sid; sid++; @@ -392,23 +508,22 @@ nconformelements = (int)this->sid2index.size(); nconformvertices = (int)sid; - newelements = xNew(nconformelements*this->elementswidth); + newelements = xNew(nconformelements*this->GetNumberOfNodes()); newmeshX = xNew(nconformvertices); newmeshY = xNew(nconformvertices); - for(int i=0;i coords(3,0.); - gmesh->NodeVec()[i].GetCoordinates(coords); - newmeshX[sid] = coords[0]; - newmeshY[sid] = coords[1]; - } + for(int i=0;i coords(3,0.); + gmesh->NodeVec()[i].GetCoordinates(coords); + newmeshX[sid] = coords[0]; + newmeshY[sid] = coords[1]; } for(int i=0;isid2index.size();i++){//over the sid (fill the ISSM elements) - for(int j=0;jelementswidth;j++) { + for(int j=0;jGetNumberOfNodes();j++) { geoel = gmesh->ElementVec()[this->sid2index[i]]; sid = vertex_index2sid[geoel->NodeIndex(j)]; - newelements[i*this->elementswidth+j]=(int)sid+1;//C to Matlab indexing + newelements[i*this->GetNumberOfNodes()+j]=(int)sid+1;//C to Matlab indexing } /*Verify the Jacobian determinant. If detJ<0, swap the 2 first postions: @@ -418,7 +533,7 @@ 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; + a=newelements[i*this->GetNumberOfNodes()+0]-1; + b=newelements[i*this->GetNumberOfNodes()+1]-1; + c=newelements[i*this->GetNumberOfNodes()+2]-1; xa=newmeshX[a]; ya=newmeshY[a]; @@ -430,12 +545,12 @@ /*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 + newelements[i*this->GetNumberOfNodes()+0]=b+1;//a->b + newelements[i*this->GetNumberOfNodes()+1]=a+1;//b->a } } /*Setting outputs*/ - nvertices = nconformvertices; - nelements = nconformelements; + *nvertices = nconformvertices; + *nelements = nconformelements; *px = newmeshX; *py = newmeshY; @@ -446,54 +561,5 @@ } /*}}}*/ -void AdaptiveMeshRefinement::FindElements(int &numberofpoints,double* xp,double* yp,TPZGeoMesh *gmesh,int &hlevel,std::vector &elements){/*{{{*/ - - if(!gmesh) _error_("Impossible to set elements: gmesh is NULL!\n"); - - if(false){ - this->AllElements(gmesh,elements); //uniform, refine all elements! - return; - } - - /*Intermediaries*/ - elements.clear(); - double D1 = this->regionlevel1; - double Dhmax = this->regionlevelmax; - int hmax = this->levelmax; - double alpha = (hmax==1) ? 0. : log(D1/Dhmax)/(hmax-1.); - double Di = D1/exp(alpha*(hlevel-1)); - int side2D = 6; - double distance,value; - - /*Find elements near the points */ - for(int i=0;iNElements();i++){ - if(gmesh->Element(i)->MaterialId()!=this->GetElemMaterialID()) continue; - if(gmesh->Element(i)->HasSubElement()) continue; - if(gmesh->Element(i)->Level()>=hlevel) continue; - TPZVec qsi(2,0.); - TPZVec centerPoint(3,0.); - gmesh->Element(i)->CenterPoint(side2D, qsi); - gmesh->Element(i)->X(qsi, centerPoint); - distance = Di; - for (int j=0;j &elements){/*{{{*/ - /* Uniform refinement. This refines the entire mesh */ - int nelements = gmesh->NElements(); - elements.clear(); - for(int i=0;iElement(i)->MaterialId()!=this->GetElemMaterialID()) continue; - if(gmesh->Element(i)->HasSubElement()) continue; - elements.push_back(i); - } -} -/*}}}*/ -void AdaptiveMeshRefinement::CreateInitialMesh(int &nvertices,int &nelements,int &width,double* x,double* y,int* elements){/*{{{*/ +void AdaptiveMeshRefinement::CreateInitialMesh(int &nvertices,int &nelements,double* x,double* y,int* elements){/*{{{*/ /* IMPORTANT! elements come in Matlab indexing @@ -502,8 +568,7 @@ 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"); - this->SetElementWidth(width); /*Verify and creating initial mesh*/ - if(this->fathermesh || this->currentmesh) _error_("Initial mesh already exists!"); + if(this->fathermesh || this->previousmesh) _error_("Initial mesh already exists!"); this->fathermesh = new TPZGeoMesh(); @@ -525,12 +590,12 @@ long index; const int mat = this->GetElemMaterialID(); - TPZManVector elem(this->elementswidth,0); + TPZManVector elem(this->GetNumberOfNodes(),0); this->sid2index.clear(); for(int i=0;ielementswidth;j++) elem[j]=elements[i*this->elementswidth+j]-1;//Convert Matlab to C indexing + for(int j=0;jGetNumberOfNodes();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->elementswidth){ + switch(this->GetNumberOfNodes()){ case 3: this->fathermesh->CreateGeoElement(ETriangle,elem,mat,index,reftype); break; default: _error_("mesh not supported yet"); @@ -543,6 +608,6 @@ /*Build element and node connectivities*/ this->fathermesh->BuildConnectivity(); - /*Set current mesh*/ - this->currentmesh=new TPZGeoMesh(*this->fathermesh); + /*Set previous mesh*/ + this->previousmesh=new TPZGeoMesh(*this->fathermesh); } /*}}}*/ @@ -629,24 +694,9 @@ } /*}}}*/ -void AdaptiveMeshRefinement::SetLevelMax(int &h){/*{{{*/ - this->levelmax = h; -} -/*}}}*/ -void AdaptiveMeshRefinement::SetRegions(double &D1,double Dhmax){/*{{{*/ - this->regionlevel1 = D1; - this->regionlevelmax = Dhmax; -} -/*}}}*/ -void AdaptiveMeshRefinement::SetElementWidth(int &width){/*{{{*/ - if(width!=3) _error_("elementswidth not supported yet!"); - this->elementswidth = width; -} -/*}}}*/ -void AdaptiveMeshRefinement::CheckMesh(int &nvertices,int &nelements,int &width,double** px,double** py,int** pelements){/*{{{*/ +void AdaptiveMeshRefinement::CheckMesh(int* nvertices,int* nelements,double** px,double** py,int** pelements){/*{{{*/ /*Basic verification*/ if(nvertices<=0) _error_("Impossible to continue: nvertices <=0!\n"); if(nelements<=0) _error_("Impossible to continue: nelements <=0!\n"); - if(width!=3) _error_("Impossible to continue: width !=3!\n"); if(!px) _error_("Impossible to continue: px is NULL!\n"); if(!py) _error_("Impossible to continue: py is NULL!\n"); @@ -656,19 +706,20 @@ std::set elemvertices; elemvertices.clear(); - for(int i=0;iGetNumberOfNodes();j++) { + elemvertices.insert((*pelements)[i*this->GetNumberOfNodes()+j]); + } + } + if(elemvertices.size()!=*nvertices) _error_("Impossible to continue: elemvertices.size() != nvertices!\n"); //Verify if there are inf or NaN in coords - for(int i=0;iGetNumberOfNodes();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"); } } @@ -676,2 +727,152 @@ } /*}}}*/ +void AdaptiveMeshRefinement::PrintGMeshVTK(TPZGeoMesh* gmesh,std::ofstream &file,bool matColor){/*{{{*/ + + file.clear(); + long nelements = gmesh->NElements(); + TPZGeoEl *gel; + std::stringstream node, connectivity, type, material; + + //Header + file << "# vtk DataFile Version 3.0" << std::endl; + file << "TPZGeoMesh VTK Visualization" << std::endl; + file << "ASCII" << std::endl << std::endl; + file << "DATASET UNSTRUCTURED_GRID" << std::endl; + file << "POINTS "; + + long actualNode = -1, size = 0, nVALIDelements = 0; + for(long el = 0; el < nelements; el++){ + gel = gmesh->ElementVec()[el]; + if(!gel )//|| (gel->Type() == EOned && !gel->IsLinearMapping()))//Exclude Arc3D and Ellipse3D + { + continue; + } + if(gel->HasSubElement()) + { + continue; + } + MElementType elt = gel->Type(); + int elNnodes = MElementType_NNodes(elt); + + size += (1+elNnodes); + connectivity << elNnodes; + + for(int t = 0; t < elNnodes; t++) + { + for(int c = 0; c < 3; c++) + { + double coord = gmesh->NodeVec()[gel->NodeIndex(t)].Coord(c); + node << coord << " "; + } + node << std::endl; + + actualNode++; + connectivity << " " << actualNode; + } + connectivity << std::endl; + + int elType = this->GetVTK_ElType(gel); + type << elType << std::endl; + + if(matColor == true) + { + material << gel->MaterialId() << std::endl; + } + else + { + material << gel->Index() << std::endl; + } + + nVALIDelements++; + } + node << std::endl; + actualNode++; + file << actualNode << " float" << std::endl << node.str(); + + file << "CELLS " << nVALIDelements << " "; + + file << size << std::endl; + file << connectivity.str() << std::endl; + + file << "CELL_TYPES " << nVALIDelements << std::endl; + file << type.str() << std::endl; + + file << "CELL_DATA" << " " << nVALIDelements << std::endl; + file << "FIELD FieldData 1" << std::endl; + if(matColor == true) + { + file << "material 1 " << nVALIDelements << " int" << std::endl; + } + else + { + file << "ElementIndex 1 " << nVALIDelements << " int" << std::endl; + } + file << material.str(); + file.close(); +} +/*}}}*/ +int AdaptiveMeshRefinement::GetVTK_ElType(TPZGeoEl * gel){/*{{{*/ + + MElementType pzElType = gel->Type(); + + int elType = -1; + switch (pzElType) + { + case(EPoint): + { + elType = 1; + break; + } + case(EOned): + { + elType = 3; + break; + } + case (ETriangle): + { + elType = 5; + break; + } + case (EQuadrilateral): + { + elType = 9; + break; + } + case (ETetraedro): + { + elType = 10; + break; + } + case (EPiramide): + { + elType = 14; + break; + } + case (EPrisma): + { + elType = 13; + break; + } + case (ECube): + { + elType = 12; + break; + } + default: + { + std::cout << "Element type not found on " << __PRETTY_FUNCTION__ << std::endl; + DebugStop(); + break; + } + } + if(elType == -1) + { + std::cout << "Element type not found on " << __PRETTY_FUNCTION__ << std::endl; + std::cout << "MIGHT BE CURVED ELEMENT (quadratic or quarter point)" << std::endl; + DebugStop(); + } + + return elType; +} +/*}}}*/ + Index: /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h =================================================================== --- /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h (revision 22099) +++ /issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h (revision 22100) @@ -37,4 +37,5 @@ #include #include +#include #include @@ -46,6 +47,22 @@ public: - - /*Public methods*/ + /*Public attributes{{{*/ + /* Filter (distance to the points) + * to refine: + * radius_h = initial_radius * gradation ^ (level_max-h) + * to unirefine + * radius_h = lag * initial_radius * gradation ^ (level_max-h) + */ + int level_max; //max level of refinement + double radius_level_max; //initial radius which in the elements will be refined with level max + double gradation; //geometric progression ratio to calculate radius of level h + double lag; //lag used in the unrefine process + /*Target and estimators*/ + double groundingline_distance; //if groundingline_distance>initial_radius, groundingline_distance will be used instead initial_radius + double icefront_distance; //if icefront_distance>initial_radius, icefront_distance will be used instead initial_radius + double thicknesserror_threshold; //if ==0, it will not be used + double deviatoricerror_threshold; //if ==0, it will not be used + /*}}}*/ + /*Public methods{{{*/ /* Constructor, destructor etc*/ AdaptiveMeshRefinement(); @@ -53,40 +70,33 @@ AdaptiveMeshRefinement & operator= (const AdaptiveMeshRefinement &cp); virtual ~AdaptiveMeshRefinement(); - /*General methods*/ void CleanUp(); void Initialize(); - void SetLevelMax(int &h); - void SetRegions(double &D1,double Dhmax); - void SetElementWidth(int &width); - void Execute(bool &amr_verbose,int &numberofelements, - double* partiallyfloatedelements,double *masklevelset,double* deviatorictensorerror,double* thicknesserror, - int &newnumberofvertices,int &newnumberofelements,double** x,double** y,int** elementslist); - void CreateInitialMesh(int &nvertices,int &nelements,int &width,double* x,double* y,int* elements); - TPZGeoMesh* CreateRefPatternMesh(TPZGeoMesh* gmesh); - void CheckMesh(int &nvertices,int &nelements,int &width,double** px,double** py,int** pelements); - + 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 CreateInitialMesh(int &nvertices,int &nelements,double* x,double* y,int* elements); + void CheckMesh(int* nvertices,int* nelements,double** px,double** py,int** pelements); + /*}}}*/ private: - /*Private attributes*/ - 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 + /*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 *currentmesh; // Current Mesh is the refined mesh - - /*Private methods*/ - void RefinementProcess(bool &amr_verbose,double* partiallyfloatedelements,double* masklevelset,double* deviatorictensorerror,double* thicknesserror); + TPZGeoMesh *previousmesh; // Previous Mesh is the refined 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(TPZGeoMesh *gmesh); - void GetMesh(int &nvertices,int &nelements,double** px,double** py,int** pelements); - void FindElements(int &numberofpoints,double* xp,double* yp,TPZGeoMesh *gmesh,int &hlevel,std::vector &elements); - void AllElements(TPZGeoMesh *gmesh,std::vector &elements); - inline int GetElemMaterialID(){return 1;} + 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); + TPZGeoMesh* CreateRefPatternMesh(TPZGeoMesh* gmesh); + inline int GetElemMaterialID(){return 1;} + inline int GetNumberOfNodes(){return 3;} + void PrintGMeshVTK(TPZGeoMesh *gmesh,std::ofstream &file,bool matColor=true); + int GetVTK_ElType(TPZGeoEl* gel); + /*}}}*/ }; Index: /issm/trunk-jpl/src/c/classes/AmrBamg.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/AmrBamg.cpp (revision 22099) +++ /issm/trunk-jpl/src/c/classes/AmrBamg.cpp (revision 22100) @@ -137,5 +137,5 @@ /*verify if the metric will be reseted or not*/ - if(!this->keepmetric){ + if(this->keepmetric==0){ if(this->options->metric) xDelete(this->options->metric); this->options->metricSize[0] = 0; Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 22099) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 22100) @@ -3590,5 +3590,4 @@ xDelete(elementlevelset); delete vmasklevelset; - } /*}}}*/ @@ -3632,4 +3631,90 @@ delete vxc; delete vyc; +} +/*}}}*/ +void FemModel::GetZeroLevelSetPoints(IssmDouble** pzerolevelset_points,int &numberofpoints,int levelset_type){/*{{{*/ + + /*Here, "zero level set" means grounding line or ice front, depending on the level set type*/ + /*pzerolevelset_points are the element center points with zero level set. X and Y coords*/ + if(levelset_type!=MaskGroundediceLevelsetEnum && levelset_type!=MaskIceLevelsetEnum){ + _error_("level set type not implemented yet!"); + } + + /*Outputs*/ + IssmDouble* zerolevelset_points = NULL; + int npoints = 0; + + /*Intermediaries*/ + int elementswidth = this->GetElementsWidth(); + int numberofelements = this->elements->NumberOfElements(); + int* elem_vertices = xNew(elementswidth); + IssmDouble* levelset = xNew(elementswidth); + 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; + IssmDouble xcenter,ycenter; + + /*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(); + + /*Find the number of points*/ + npoints=0; + for(int i=0;i(x_zerolevelset[i])) npoints++; + + /*Keep just the element center coordinates with zero level set (compact the structure)*/ + zerolevelset_points=xNew(2*npoints);//x and y + count=0; + for(int i=0;i(x_zerolevelset[i])){ + zerolevelset_points[2*count] = x_zerolevelset[i]; + zerolevelset_points[2*count+1] = y_zerolevelset[i]; + count++; + } + } + + /*Assign outputs*/ + numberofpoints = npoints; + (*pzerolevelset_points) = zerolevelset_points; + + /*Cleanup*/ + xDelete(elem_vertices); + xDelete(levelset); + xDelete(x_zerolevelset); + xDelete(y_zerolevelset); + xDelete(x); + xDelete(y); + xDelete(z); + delete vx_zerolevelset; + delete vy_zerolevelset; } /*}}}*/ @@ -4586,24 +4671,219 @@ /*Intermediaries*/ + int numberofvertices = this->vertices->NumberOfVertices(); + IssmDouble* levelset_points= NULL; + IssmDouble* x = NULL; + IssmDouble* y = NULL; + IssmDouble* z = NULL; + int numberofpoints; + IssmDouble distance; + + /*Get vertices coordinates*/ + VertexCoordinatesx(&x,&y,&z,this->vertices,false) ; + + /*Get points which level set is zero (center of elements with zero level set)*/ + this->GetZeroLevelSetPoints(&levelset_points,numberofpoints,levelset_type); + + /*Find the minimal vertice distance to the zero levelset (grounding line or ice front)*/ + verticedistance=xNew(numberofvertices); + for(int i=0;i(levelset_points); + xDelete(x); + xDelete(y); + xDelete(z); +} +/*}}}*/ +#endif + +#if defined(_HAVE_NEOPZ_) && !defined(_HAVE_ADOLC_) +void FemModel::ReMeshNeopz(int* pnewnumberofvertices,int* pnewnumberofelements,IssmDouble** pnewx,IssmDouble** pnewy,IssmDouble** pnewz,int** pnewelementslist){/*{{{*/ + + /*pnewelementslist keep vertices in Matlab indexing*/ + int my_rank = IssmComm::GetRank(); + int numberofelements = this->elements->NumberOfElements(); + IssmDouble* element_label = xNewZeroInit(numberofelements); + int numberofpoints = -1; + IssmDouble* xylist = NULL; + IssmDouble* newx = NULL; + IssmDouble* newy = NULL; + IssmDouble* newz = NULL; + int* newelementslist = NULL; + int newnumberofvertices = -1; + 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); + + if(my_rank==0){ + this->amr->ExecuteRefinement(numberofpoints,xylist,&newnumberofvertices,&newnumberofelements,&newx,&newy,&newelementslist); + newz=xNewZeroInit(newnumberofvertices); + if(newnumberofvertices<=0 || newnumberofelements<=0) _error_("Error in the ReMeshNeopz."); + } + else{ + newx=xNew(newnumberofvertices); + newy=xNew(newnumberofvertices); + newz=xNew(newnumberofvertices); + newelementslist=xNew(newnumberofelements*this->GetElementsWidth()); + } + + /*Send new mesh to others CPU*/ + ISSM_MPI_Bcast(&newnumberofvertices,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(&newnumberofelements,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm()); + + /*Assign the pointers*/ + (*pnewelementslist) = newelementslist; //Matlab indexing + (*pnewx) = newx; + (*pnewy) = newy; + (*pnewz) = newz; + *pnewnumberofvertices= newnumberofvertices; + *pnewnumberofelements= newnumberofelements; + + /*Cleanup*/ + xDelete(element_label); + xDelete(xylist); +} +/*}}}*/ +void FemModel::InitializeAdaptiveRefinementNeopz(void){/*{{{*/ + + /*Define variables*/ + int my_rank = IssmComm::GetRank(); + int numberofvertices = this->vertices->NumberOfVertices(); + int numberofelements = this->elements->NumberOfElements(); + IssmDouble* x = NULL; + IssmDouble* y = NULL; + IssmDouble* z = NULL; + int* elements = NULL; + int level_max = -1; + IssmDouble radius_level_max = -1; + IssmDouble gradation = -1; + IssmDouble lag = -1; + IssmDouble groundingline_distance = -1; + IssmDouble icefront_distance = -1; + IssmDouble thicknesserror_threshold = -1; + IssmDouble deviatoricerror_threshold = -1; + + /*Initialize field as NULL for now*/ + this->amr = NULL; + + /*Get vertices coordinates of the coarse mesh (father mesh)*/ + /*elements comes in Matlab indexing*/ + this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elements); + + /*Get amr parameters*/ + this->parameters->FindParam(&level_max,AmrLevelMaxEnum); + this->parameters->FindParam(&radius_level_max,AmrRadiusLevelMaxEnum); + this->parameters->FindParam(&gradation,AmrGradationEnum); + this->parameters->FindParam(&lag,AmrLagEnum); + this->parameters->FindParam(&groundingline_distance,AmrGroundingLineDistanceEnum); + this->parameters->FindParam(&icefront_distance,AmrIceFrontDistanceEnum); + this->parameters->FindParam(&thicknesserror_threshold,AmrThicknessErrorThresholdEnum); + this->parameters->FindParam(&deviatoricerror_threshold,AmrDeviatoricErrorThresholdEnum); + + /*Create initial mesh (coarse mesh) in neopz data structure*/ + /*Just CPU #0 should keep AMR object*/ + /*Initialize refinement pattern*/ + this->SetRefPatterns(); + this->amr = new AdaptiveMeshRefinement(); + this->amr->level_max = level_max; + this->amr->radius_level_max = radius_level_max; + this->amr->gradation = gradation; + this->amr->lag = lag; + this->amr->groundingline_distance = groundingline_distance; + this->amr->icefront_distance = icefront_distance; + this->amr->thicknesserror_threshold = thicknesserror_threshold; + this->amr->deviatoricerror_threshold = deviatoricerror_threshold; + if(my_rank==0){ + this->amr->CreateInitialMesh(numberofvertices,numberofelements,x,y,elements); + } + + /*Free the vectors*/ + xDelete(x); + xDelete(y); + xDelete(z); + xDelete(elements); +} +/*}}}*/ +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 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) ; - + 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++){ @@ -4611,171 +4891,112 @@ element->GetInputListOnVertices(levelset,levelset_type); element->GetVerticesSidList(elem_vertices); - sid = element->Sid(); - xcenter = NAN; - ycenter = NAN; + 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,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); + velement_label->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){/*{{{*/ + + /*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* elementdistance; + + /*Intermediaries*/ + int numberofelements = this->elements->NumberOfElements(); + IssmDouble* levelset_points= NULL; + IssmDouble* xc = NULL; + IssmDouble* yc = NULL; + 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)*/ + 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;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; -} -/*}}}*/ -#endif - -#if defined(_HAVE_NEOPZ_) && !defined(_HAVE_ADOLC_) -void FemModel::ReMeshNeopz(int* pnewnumberofvertices,int* pnewnumberofelements,IssmDouble** pnewx,IssmDouble** pnewy,IssmDouble** pnewz,int** pnewelementslist){/*{{{*/ - - /*pnewelementslist keep vertices in Matlab indexing*/ - int my_rank = IssmComm::GetRank(); - bool amr_verbose = true; //itapopo - IssmDouble* x = NULL; - IssmDouble* y = NULL; - IssmDouble* z = NULL; - int* elementslist = NULL; - int oldnumberofelements = this->elements->NumberOfElements(); - int newnumberofvertices = -1; - int newnumberofelements = -1; - IssmDouble* partiallyfloatedelements= NULL;//itapopo verify if it will be used - IssmDouble* masklevelset = NULL; - IssmDouble* deviatorictensorerror = NULL; - IssmDouble* thicknesserror = NULL; - - /*Get the elements in which grounding line goes through*/ - //itapopo verificar se irá usar esse this->GetPartiallyFloatedElements(&partiallyfloatedelements); - /*Get mean mask level set over each element*/ - this->MeanGroundedIceLevelSet(&masklevelset); - /*Get the deviatoric error estimator*/ - this->ZZErrorEstimator(&deviatorictensorerror); - /*Get the thickness error estimator*/ - this->ThicknessZZErrorEstimator(&thicknesserror); - - if(my_rank==0){ - this->amr->Execute(amr_verbose,oldnumberofelements,partiallyfloatedelements,masklevelset,deviatorictensorerror,thicknesserror, - newnumberofvertices,newnumberofelements,&x,&y,&elementslist); - z=xNewZeroInit(newnumberofvertices); - if(newnumberofvertices<=0 || newnumberofelements<=0) _error_("Error in the ReMeshNeopz."); - } - else{ - x=xNew(newnumberofvertices); - y=xNew(newnumberofvertices); - z=xNew(newnumberofvertices); - elementslist=xNew(newnumberofelements*this->GetElementsWidth()); - } - - /*Send new mesh to others CPU*/ - ISSM_MPI_Bcast(&newnumberofvertices,1,ISSM_MPI_INT,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(&newnumberofelements,1,ISSM_MPI_INT,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(x,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(y,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(z,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(elementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm()); - - /*Assign the pointers*/ - (*pnewelementslist) = elementslist; //Matlab indexing - (*pnewx) = x; - (*pnewy) = y; - (*pnewz) = z; - *pnewnumberofvertices= newnumberofvertices; - *pnewnumberofelements= newnumberofelements; - - /*Cleanup*/ - xDelete(masklevelset); - xDelete(deviatorictensorerror); - xDelete(thicknesserror); -} -/*}}}*/ -void FemModel::InitializeAdaptiveRefinementNeopz(void){/*{{{*/ - - /*Define variables*/ - int my_rank = IssmComm::GetRank(); - int numberofvertices = this->vertices->NumberOfVertices(); - int numberofelements = this->elements->NumberOfElements(); - IssmDouble* x = NULL; - IssmDouble* y = NULL; - IssmDouble* z = NULL; - int* elements = NULL; - int elementswidth = this->GetElementsWidth(); //just tria elements in this version. Itapopo: - int levelmax = 0; - IssmDouble regionlevel1 = 0.; - IssmDouble regionlevelmax = 0.; - - /*Initialize field as NULL for now*/ - this->amr = NULL; - - /*Get vertices coordinates of the coarse mesh (father mesh)*/ - /*elements comes in Matlab indexing*/ - this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elements); - - /*Get amr parameters*/ - this->parameters->FindParam(&levelmax,AmrLevelMaxEnum); - this->parameters->FindParam(®ionlevel1,AmrRegionLevel1Enum); - this->parameters->FindParam(®ionlevelmax,AmrRegionLevelMaxEnum); - - /*Create initial mesh (coarse mesh) in neopz data structure*/ - /*Just CPU #0 should keep AMR object*/ - /*Initialize refinement pattern*/ - this->SetRefPatterns(); - if(my_rank==0){ - 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); - } - - /*Free the vectors*/ - xDelete(x); - xDelete(y); - xDelete(z); - xDelete(elements); + xDelete(levelset_points); + xDelete(xc); + xDelete(yc); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 22099) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 22100) @@ -178,4 +178,5 @@ void MeanGroundedIceLevelSet(IssmDouble** pmasklevelset); void GetElementCenterCoordinates(IssmDouble** pxc,IssmDouble** pyc); + void GetZeroLevelSetPoints(IssmDouble** pzerolevelset_points,int &numberofpoints,int levelset_type); #endif @@ -191,4 +192,8 @@ 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); #endif Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 22099) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 22100) @@ -128,6 +128,11 @@ case AmrNeopzEnum: parameters->AddObject(iomodel->CopyConstantObject("md.amr.level_max",AmrLevelMaxEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.amr.region_level_1",AmrRegionLevel1Enum)); - parameters->AddObject(iomodel->CopyConstantObject("md.amr.region_level_max",AmrRegionLevelMaxEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.radius_level_max",AmrRadiusLevelMaxEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.gradation",AmrGradationEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.lag",AmrLagEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.groundingline_distance",AmrGroundingLineDistanceEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.icefront_distance",AmrIceFrontDistanceEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.thicknesserror_threshold",AmrThicknessErrorThresholdEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.amr.deviatoricerror_threshold",AmrDeviatoricErrorThresholdEnum)); break; #endif Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 22099) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 22100) @@ -861,6 +861,6 @@ AmrNeopzEnum, AmrLevelMaxEnum, - AmrRegionLevel1Enum, - AmrRegionLevelMaxEnum, + AmrLagEnum, + AmrRadiusLevelMaxEnum, AmrBamgEnum, AmrHminEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 22099) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 22100) @@ -834,6 +834,6 @@ case AmrNeopzEnum : return "AmrNeopz"; case AmrLevelMaxEnum : return "AmrLevelMax"; - case AmrRegionLevel1Enum : return "AmrRegionLevel1"; - case AmrRegionLevelMaxEnum : return "AmrRegionLevelMax"; + case AmrLagEnum : return "AmrLag"; + case AmrRadiusLevelMaxEnum : return "AmrRadiusLevelMax"; case AmrBamgEnum : return "AmrBamg"; case AmrHminEnum : return "AmrHmin"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 22099) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 22100) @@ -852,6 +852,6 @@ else if (strcmp(name,"AmrNeopz")==0) return AmrNeopzEnum; else if (strcmp(name,"AmrLevelMax")==0) return AmrLevelMaxEnum; - else if (strcmp(name,"AmrRegionLevel1")==0) return AmrRegionLevel1Enum; - else if (strcmp(name,"AmrRegionLevelMax")==0) return AmrRegionLevelMaxEnum; + else if (strcmp(name,"AmrLag")==0) return AmrLagEnum; + else if (strcmp(name,"AmrRadiusLevelMax")==0) return AmrRadiusLevelMaxEnum; else if (strcmp(name,"AmrBamg")==0) return AmrBamgEnum; else if (strcmp(name,"AmrHmin")==0) return AmrHminEnum; Index: /issm/trunk-jpl/src/m/classes/amr.js =================================================================== --- /issm/trunk-jpl/src/m/classes/amr.js (revision 22099) +++ /issm/trunk-jpl/src/m/classes/amr.js (revision 22100) @@ -7,20 +7,34 @@ //methods this.setdefaultparameters = function(){// {{{ - //level_max: 2 to 4 - this.level_max=2; - - //region_level_1: region around (m) the discontinuity (grounding line or ice front) where the mesh will be refined once (h=1). - this.region_level_1=20000.; - - //region_level_max: region around (m) the discontinuity (grounding line or ice front) where the mesh will be refined with max level of refinement (h=level_max). - this.region_level_max=15000.; + this.hmin = 100.; + this.hmax = 100.e3; + this.fieldname = "Vel"; + this.err = 3.; + this.keepmetric = 1; + this.gradation = 1.5; + this.groundingline_resolution = 500.; + this.groundingline_distance = 0; + this.icefront_resolution = 500; + this.icefront_distance = 0; + this.thicknesserror_resolution = 500; + this.thicknesserror_threshold = 0; + this.deviatoricerror_resolution = 500; + this.deviatoricerror_threshold = 0; }// }}} this.disp= function(){// {{{ - console.log(sprintf(' amr parameters:')); - fielddisplay(this,'level_max','maximum refinement level (1, 2, 3 or 4)'); - fielddisplay(this,'region_level_1','region which will be refined once (level 1) [ m ]'); - fielddisplay(this,'region_level_max','region which will be refined with level_max [ m ]'); - + fielddisplay(this,'hmin','minimum element length'); + fielddisplay(this,'hmax','maximum element length'); + fielddisplay(this,'fieldname','name of input that will be used to compute the metric (should be an input of FemModel)'); + fielddisplay(this,'keepmetric','indicates whether the metric should be kept every remeshing time'); + fielddisplay(this,'gradation','maximum ratio between two adjacent edges'); + fielddisplay(this,'groundingline_resolution','element length near the grounding line'); + fielddisplay(this,'groundingline_distance','distance around the grounding line which elements will be refined'); + fielddisplay(this,'icefront_resolution','element length near the ice front'); + fielddisplay(this,'icefront_distance','distance around the ice front which elements will be refined'); + fielddisplay(this,'thicknesserror_resolution','element length when thickness error estimator is used'); + fielddisplay(this,'thicknesserror_threshold','maximum threshold 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'); }// }}} this.classname= function(){// {{{ @@ -29,18 +43,33 @@ }// }}} this.checkconsistency = function(md,solution,analyses) { //{{{ - - checkfield(md,'fieldname','amr.level_max','numel',[1],'>=',0,'<=',4); - checkfield(md,'fieldname','amr.region_level_1','numel',[1],'>',0,'NaN',1,'Inf',1); - checkfield(md,'fieldname','amr.region_level_max','numel',[1],'>',0,'NaN',1,'Inf',1); - if (this.region_level_1-this.region_level_max<0.2*this.region_level_1){ - md.checkmessage('region_level_max should be lower than 80% of region_level_1'); - } + checkfield(md,'fieldname','amr.hmax','numel',[1],'>',0,'NaN',1); + checkfield(md,'fieldname','amr.hmin','numel',[1],'>',0,'<',this.hmax,'NaN',1); + checkfield(md,'fieldname','amr.keepmetric','numel',[1],'>=',0,'<=',1,'NaN',1); + checkfield(md,'fieldname','amr.gradation','numel',[1],'>=',1.1,'<=',5,'NaN',1); + checkfield(md,'fieldname','amr.groundingline_resolution','numel',[1],'>',0,'<',this.hmax,'NaN',1); + checkfield(md,'fieldname','amr.groundingline_distance','numel',[1],'>=',0,'NaN',1,'Inf',1); + checkfield(md,'fieldname','amr.icefront_resolution','numel',[1],'>',0,'<',this.hmax,'NaN',1); + checkfield(md,'fieldname','amr.icefront_distance','numel',[1],'>=',0,'NaN',1,'Inf',1); + 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.deviatoricerror_resolution','numel',[1],'>',0,'<',this.hmax,'NaN',1); + checkfield(md,'fieldname','amr.deviatoricerror_threshold','numel',[1],'>=',0,'<=',1,'NaN',1); } // }}} this.marshall=function(md,prefix,fid) { //{{{ - - WriteData(fid,prefix,'object',this,'fieldname','level_max','format','Integer'); - WriteData(fid,prefix,'object',this,'fieldname','region_level_1','format','Double'); - WriteData(fid,prefix,'object',this,'fieldname','region_level_max','format','Double'); - + WriteData(fid,prefix,'name','md.amr.type','data',1,'format','Integer'); + WriteData(fid,prefix,'object',this,'fieldname','hmin','format','Double'); + WriteData(fid,prefix,'object',this,'fieldname','hmax','format','Double'); + WriteData(fid,prefix,'object',this,'fieldname','fieldname','format','String'); + WriteData(fid,prefix,'object',this,'fieldname','err','format','Double'); + WriteData(fid,prefix,'object',this,'fieldname','keepmetric','format','Integer'); + WriteData(fid,prefix,'object',this,'fieldname','gradation','format','Double'); + WriteData(fid,prefix,'object',this,'fieldname','groundingline_resolution','format','Double'); + WriteData(fid,prefix,'object',this,'fieldname','groundingline_distance','format','Double'); + WriteData(fid,prefix,'object',this,'fieldname','icefront_resolution','format','Double'); + WriteData(fid,prefix,'object',this,'fieldname','icefront_distance','format','Double'); + 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','deviatoricerror_resolution','format','Double'); + WriteData(fid,prefix,'object',this,'fieldname','deviatoricerror_threshold','format','Double'); }//}}} this.fix=function() { //{{{ @@ -48,7 +77,18 @@ //properties // {{{ - this.level_max = 0; - this.region_level_1 = 0.; - this.region_level_max = 0.; + this.hmin = 0.; + this.hmax = 0.; + this.fieldname = ""; + this.err = 0.; + this.keepmetric = 0; + this.gradation = 0.; + this.groundingline_resolution = 0.; + this.groundingline_distance = 0.; + this.icefront_resolution = 0.; + this.icefront_distance = 0.; + this.thicknesserror_resolution = 0.; + this.thicknesserror_threshold = 0.; + this.deviatoricerror_resolution = 0.; + this.deviatoricerror_threshold = 0.; this.setdefaultparameters(); Index: /issm/trunk-jpl/src/m/classes/amr.m =================================================================== --- /issm/trunk-jpl/src/m/classes/amr.m (revision 22099) +++ /issm/trunk-jpl/src/m/classes/amr.m (revision 22100) @@ -2,11 +2,55 @@ % % Usage: -% amr=amr(); +% md.amr=amr(); classdef amr properties (SetAccess=public) - level_max = 0; - region_level_1 = 0; - region_level_max = 0; + hmin = 0.; + hmax = 0.; + fieldname = ''; + err = 0.; + keepmetric = 0; + gradation = 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 (Static) + function self = loadobj(self) % {{{ + % This function is directly called by matlab when a model object is + % loaded. Update old properties here + + if verLessThan('matlab','7.9'), + disp('Warning: your matlab version is old and there is a risk that load does not work correctly'); + disp(' if the model is not loaded correctly, rename temporarily loadobj so that matlab does not use it'); + + % This is a Matlab bug: all the fields of md have their default value + % Example of error message: + % Warning: Error loading an object of class 'model': + % Undefined function or method 'exist' for input arguments of type 'cell' + % + % This has been fixed in MATLAB 7.9 (R2009b) and later versions + end + + %2017 September 15th + if isstruct(self), + disp('WARNING: updating amr. Now the default is amr with bamg'); + disp(' some old fields were not converted'); + disp(' see the new fields typing md.amr and modify them properly'); + obj2 = self; + self = amr(); + %Converting region_level_max to groundingline_distance + if(obj2.region_level_max>0 && obj2.region_level_1>obj2.region_level_max) + self.groundingline_distance = obj2.region_level_max; + self.keepmetric = 0; + self.fieldname = 'None'; + end + end + end% }}} end methods @@ -21,43 +65,80 @@ function self = setdefaultparameters(self) % {{{ - %level_max: 2 to 4 - self.level_max=2; + %hmin and hmax + self.hmin=100.; + self.hmax=100.e3; - %region_level_1: region around (m) the discontinuity (grounding line or ice front) where the mesh will be refined once (h=1). - self.region_level_1=20000.; - - %region_level_max: region around (m) the discontinuity (grounding line or ice front) where the mesh will be refined with max level of refinement (h=level_max). - self.region_level_max=15000.; + %fields + self.fieldname ='Vel'; + self.err=3.; + + %keep metric? + self.keepmetric=1; + + %control of element lengths + self.gradation=1.5; + + %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) % {{{ - md = checkfield(md,'fieldname','amr.level_max','numel',[1],'>=',0,'<=',4); - md = checkfield(md,'fieldname','amr.region_level_1','numel',[1],'>',0,'NaN',1,'Inf',1); - md = checkfield(md,'fieldname','amr.region_level_max','numel',[1],'>',0,'NaN',1,'Inf',1); - if self.region_level_1-self.region_level_max<0.2*self.region_level_1, %it was adopted 20% of the region_level_1 - md = checkmessage(md,'region_level_max should be lower than 80% of region_level_1'); - end + md = checkfield(md,'fieldname','amr.hmax','numel',[1],'>',0,'NaN',1); + 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.gradation','numel',[1],'>=',1.1,'<=',5,'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) % {{{ disp(sprintf(' amr parameters:')); - fielddisplay(self,'level_max',['maximum refinement level (1, 2, 3 or 4)']); - fielddisplay(self,'region_level_1',['region which will be refined once (level 1) [ m ]']); - fielddisplay(self,'region_level_max',['region which will be refined with level_max [ m ]']); + fielddisplay(self,'hmin',['minimum element length']); + 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,'gradation',['maximum ratio between two adjacent edges']); + 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 % }}} function marshall(self,prefix,md,fid) % {{{ - WriteData(fid,prefix,'name','md.amr.type','data',2,'format','Integer'); - WriteData(fid,prefix,'object',self,'fieldname','level_max','format','Integer'); - WriteData(fid,prefix,'object',self,'fieldname','region_level_1','format','Double'); - WriteData(fid,prefix,'object',self,'fieldname','region_level_max','format','Double'); - end % }}} - function savemodeljs(self,fid,modelname) % {{{ - - writejsdouble(fid,[modelname '.amr.level_max'],self.level_max); - writejsdouble(fid,[modelname '.amr.region_level_1'],self.region_level_1); - writejsdouble(fid,[modelname '.amr.region_level_max'],self.region_level_max); + WriteData(fid,prefix,'name','md.amr.type','data',1,'format','Integer'); + WriteData(fid,prefix,'object',self,'class','amr','fieldname','hmin','format','Double'); + 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','amr','fieldname','err','format','Double'); + WriteData(fid,prefix,'object',self,'class','amr','fieldname','keepmetric','format','Integer'); + WriteData(fid,prefix,'object',self,'class','amr','fieldname','gradation','format','Double'); + 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/classes/amr.py =================================================================== --- /issm/trunk-jpl/src/m/classes/amr.py (revision 22099) +++ /issm/trunk-jpl/src/m/classes/amr.py (revision 22100) @@ -12,8 +12,18 @@ def __init__(self): # {{{ - self.level_max = 0. - self.region_level_1 = 0. - self.region_level_max = 0. - + self.hmin = 0. + self.hmax = 0. + self.fieldname ='' + self.err = 0. + self.keepmetric = 0. + self.gradation = 0. + self.groundingline_resolution = 0. + self.groundingline_distance = 0. + self.icefront_resolution = 0. + self.icefront_distance = 0. + self.thicknesserror_resolution = 0. + self.thicknesserror_threshold = 0. + self.deviatoricerror_resolution= 0. + self.deviatoricerror_threshold = 0. #set defaults self.setdefaultparameters() @@ -21,35 +31,66 @@ def __repr__(self): # {{{ string=" amr parameters:" - string="%s\n%s"%(string,fielddisplay(self,"level_max","maximum refinement level (1, 2, 3 or 4)")) - string="%s\n%s"%(string,fielddisplay(self,"region_level_1","region which will be refined once (level 1) [ m ]")) - string="%s\n%s"%(string,fielddisplay(self,"region_level_max","region which will be refined with level_max [ m ]")) + string="%s\n%s"%(string,fielddisplay(self,"hmin","minimum element length")) + string="%s\n%s"%(string,fielddisplay(self,"hmax","maximum element length")) + string="%s\n%s"%(string,fielddisplay(self,"fieldname","name of input that will be used to compute the metric (should be an input of FemModel)")) + string="%s\n%s"%(string,fielddisplay(self,"keepmetric","indicates whether the metric should be kept every remeshing time")) + string="%s\n%s"%(string,fielddisplay(self,"gradation","maximum ratio between two adjacent edges")) + string="%s\n%s"%(string,fielddisplay(self,"groundingline_resolution","element length near the grounding line")) + string="%s\n%s"%(string,fielddisplay(self,"groundingline_distance","distance around the grounding line which elements will be refined")) + string="%s\n%s"%(string,fielddisplay(self,"icefront_resolution","element length near the ice front")) + string="%s\n%s"%(string,fielddisplay(self,"icefront_distance","distance around the ice front which elements will be refined")) + 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,"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")) return string #}}} def setdefaultparameters(self): # {{{ - - #level_max: 2 to 4 - self.level_max=2 - - #region_level_1: region around (m) the discontinuity (grounding line or ice front) where the mesh will be refined once (h=1). - self.region_level_1=20000. - - #region_level_max: region around (m) the discontinuity (grounding line or ice front) where the mesh will be refined with max level of refinement (h=level_max). - self.region_level_max=15000. - + self.hmin = 100. + self.hmax = 100.e3 + self.fieldname = 'Vel' + self.err = 3. + self.keepmetric = 1 + self.gradation = 1.5 + 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 return self - #}}} + #}}} def checkconsistency(self,md,solution,analyses): # {{{ - md = checkfield(md,'fieldname','amr.level_max','numel',[1],'>=',0,'<=',4) - md = checkfield(md,'fieldname','amr.region_level_1','numel',[1],'>',0,'NaN',1,'Inf',1) - md = checkfield(md,'fieldname','amr.region_level_max','numel',[1],'>',0,'NaN',1,'Inf',1) - #it was adopted 20% of the region_level_1 - if self.region_level_1-self.region_level_max<0.2*self.region_level_1: - md.checkmessage("region_level_max should be lower than 80% of region_level_1") - + md = checkfield(md,'fieldname','amr.hmax','numel',[1],'>',0,'NaN',1) + md = checkfield(md,'fieldname','amr.hmin','numel',[1],'>',0,'<',self.hmax,'NaN',1) + md = checkfield(md,'fieldname','amr.keepmetric','numel',[1],'>=',0,'<=',1,'NaN',1); + md = checkfield(md,'fieldname','amr.gradation','numel',[1],'>=',1.1,'<=',5,'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); return md # }}} def marshall(self,prefix,md,fid): # {{{ - WriteData(fid,prefix,'object',self,'fieldname','level_max','format','Integer') - WriteData(fid,prefix,'object',self,'fieldname','region_level_1','format','Double') - WriteData(fid,prefix,'object',self,'fieldname','region_level_max','format','Double') + WriteData(fid,prefix,'name','md.amr.type','data',1,'format','Integer') + WriteData(fid,prefix,'object',self,'fieldname','hmin','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','hmax','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','fieldname','format','String'); + WriteData(fid,prefix,'object',self,'fieldname','err','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','keepmetric','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','gradation','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','groundingline_resolution','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','groundingline_distance','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','icefront_resolution','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','icefront_distance','format','Double'); + 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','deviatoricerror_resolution','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','deviatoricerror_threshold','format','Double'); # }}} Index: sm/trunk-jpl/src/m/classes/amrbamg.m =================================================================== --- /issm/trunk-jpl/src/m/classes/amrbamg.m (revision 22099) +++ (revision ) @@ -1,113 +1,0 @@ -%AMRBAMG Class definition -% -% Usage: -% amrbamg=amrbamg(); - -classdef amrbamg - properties (SetAccess=public) - hmin = 0.; - hmax = 0.; - fieldname = ''; - err = 0.; - keepmetric = 0; - gradation = 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 - function self = amrbamg(varargin) % {{{ - switch nargin - case 0 - self=setdefaultparameters(self); - otherwise - error('constructor not supported'); - end - end % }}} - function self = setdefaultparameters(self) % {{{ - - %hmin and hmax - self.hmin=100.; - self.hmax=100.e3; - - %fields - self.fieldname ='Vel'; - self.err=3.; - - %keep metric? - self.keepmetric=1; - - %control of element lengths - self.gradation=1.5; - - %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) % {{{ - - md = checkfield(md,'fieldname','amr.hmax','numel',[1],'>',0,'NaN',1); - 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.gradation','numel',[1],'>=',1.1,'<=',5,'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) % {{{ - disp(sprintf(' amrbamg parameters:')); - - fielddisplay(self,'hmin',['minimum element length']); - 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,'gradation',['maximum ratio between two adjacent edges']); - 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 % }}} - function marshall(self,prefix,md,fid) % {{{ - - WriteData(fid,prefix,'name','md.amr.type','data',1,'format','Integer'); - WriteData(fid,prefix,'object',self,'class','amr','fieldname','hmin','format','Double'); - 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','amr','fieldname','err','format','Double'); - WriteData(fid,prefix,'object',self,'class','amr','fieldname','keepmetric','format','Integer'); - WriteData(fid,prefix,'object',self,'class','amr','fieldname','gradation','format','Double'); - 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 % }}} - end -end Index: /issm/trunk-jpl/test/NightlyRun/test462.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test462.m (revision 22099) +++ /issm/trunk-jpl/test/NightlyRun/test462.m (revision 22100) @@ -11,5 +11,4 @@ md.transient.isgroundingline=0; %amr bamg settings, just field -md.amr=amrbamg(); md.amr.hmin=10000; md.amr.hmax=100000; Index: /issm/trunk-jpl/test/NightlyRun/test463.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test463.m (revision 22099) +++ /issm/trunk-jpl/test/NightlyRun/test463.m (revision 22100) @@ -11,5 +11,4 @@ md.transient.isgroundingline=0; %amr bamg settings, just grounding line -md.amr=amrbamg(); md.amr.hmin=10000; md.amr.hmax=100000; Index: /issm/trunk-jpl/test/NightlyRun/test464.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test464.m (revision 22099) +++ /issm/trunk-jpl/test/NightlyRun/test464.m (revision 22100) @@ -11,5 +11,4 @@ md.transient.isgroundingline=0; %amr bamg settings, just ice front -md.amr=amrbamg(); md.amr.hmin=10000; md.amr.hmax=100000; Index: /issm/trunk-jpl/test/NightlyRun/test465.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test465.m (revision 22099) +++ /issm/trunk-jpl/test/NightlyRun/test465.m (revision 22100) @@ -11,5 +11,4 @@ md.transient.isgroundingline=0; %amr bamg settings, field, grounding line and ice front -md.amr=amrbamg(); md.amr.hmin=20000; md.amr.hmax=100000;