![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 21919
- Timestamp:
- 08/09/17 10:26:20 (8 years ago)
- Location:
- issm/trunk-jpl/src
- Files:
-
- 13 edited
-
c/classes/AdaptiveMeshRefinement.cpp (modified) (14 diffs)
-
c/classes/AdaptiveMeshRefinement.h (modified) (5 diffs)
-
c/classes/AmrBamg.cpp (modified) (4 diffs)
-
c/classes/AmrBamg.h (modified) (1 diff)
-
c/classes/Elements/Element.cpp (modified) (1 diff)
-
c/classes/FemModel.cpp (modified) (13 diffs)
-
c/classes/FemModel.h (modified) (2 diffs)
-
c/modules/ModelProcessorx/CreateParameters.cpp (modified) (1 diff)
-
c/shared/Enum/EnumDefinitions.h (modified) (1 diff)
-
c/shared/Enum/EnumToStringx.cpp (modified) (1 diff)
-
c/shared/Enum/StringToEnumx.cpp (modified) (5 diffs)
-
m/classes/amrbamg.m (modified) (5 diffs)
-
m/mesh/bamg.m (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp
r21806 r21919 10 10 11 11 #include "./AdaptiveMeshRefinement.h" 12 #include "TPZVTKGeoMesh.h" 13 #include "../shared/shared.h" 14 #include "pzgeotriangle.h" 15 #include "pzreftriangle.h" 12 16 13 using namespace pzgeom; 17 14 … … 40 37 this->sid2index.resize(cp.sid2index.size()); 41 38 for(int i=0;i<cp.sid2index.size();i++) this->sid2index[i]=cp.sid2index[i]; 42 39 this->index2sid.clear(); 40 this->index2sid.resize(cp.index2sid.size()); 41 for(int i=0;i<cp.index2sid.size();i++) this->index2sid[i]=cp.index2sid[i]; 42 this->specialelementsindex.clear(); 43 this->specialelementsindex.resize(cp.specialelementsindex.size()); 44 for(int i=0;i<cp.specialelementsindex.size();i++) this->specialelementsindex[i]=cp.specialelementsindex[i]; 45 43 46 return *this; 44 45 47 } 46 48 /*}}}*/ 47 49 AdaptiveMeshRefinement::~AdaptiveMeshRefinement(){/*{{{*/ 48 49 bool ismismip = false;50 if(ismismip){//itapopo51 TPZFileStream fstr;52 std::stringstream ss;53 54 ss << this->levelmax;55 std::string AMRfile = "/home/santos/L" + ss.str() + "_amr.txt";56 57 fstr.OpenWrite(AMRfile.c_str());58 int withclassid = 1;59 this->Write(fstr,withclassid);60 }61 50 this->CleanUp(); 62 51 gRefDBase.clear(); … … 73 62 this->regionlevelmax = -1; 74 63 this->sid2index.clear(); 64 this->index2sid.clear(); 65 this->specialelementsindex.clear(); 75 66 } 76 67 /*}}}*/ … … 84 75 this->regionlevel1 = -1; 85 76 this->regionlevelmax = -1; 77 this->step = 0; 78 this->smooth_frequency = 1; 86 79 this->sid2index.clear(); 87 } 88 /*}}}*/ 89 int AdaptiveMeshRefinement::ClassId() const{/*{{{*/ 90 return 13829430; //Antartic area with ice shelves (km^2) 91 } 92 /*}}}*/ 93 void AdaptiveMeshRefinement::Read(TPZStream &buf,void *context){/*{{{*/ 94 95 try 96 { 97 /* Read the id context*/ 98 TPZSaveable::Read(buf,context); 99 /* Read class id*/ 100 int classid; 101 buf.Read(&classid,1); 102 /* Verify the class id*/ 103 if(classid!=this->ClassId()) _error_("AdaptiveMeshRefinement::Read: Error in restoring AdaptiveMeshRefinement!\n"); 104 /* Read simple attributes */ 105 buf.Read(&this->levelmax,1); 106 buf.Read(&this->elementswidth,1); 107 buf.Read(&this->regionlevel1,1); 108 buf.Read(&this->regionlevelmax,1); 109 /* Read vector attributes*/ 110 int size; 111 buf.Read(&size,1); 112 int* psid2index=xNew<int>(size); 113 buf.Read(psid2index,size); 114 this->sid2index.clear(); 115 this->sid2index.assign(psid2index,psid2index+size); 116 /* Read geometric mesh (father)*/ 117 TPZSaveable *sv1 = TPZSaveable::Restore(buf,0); 118 this->fathermesh = dynamic_cast<TPZGeoMesh*>(sv1); 119 /* Read geometric mesh (current)*/ 120 TPZSaveable *sv2 = TPZSaveable::Restore(buf,0); 121 this->currentmesh = dynamic_cast<TPZGeoMesh*>(sv2); 122 /* Cleanup*/ 123 xDelete<int>(psid2index); 124 } 125 catch(const std::exception& e) 126 { 127 _error_("AdaptiveMeshRefinement::Read: Exception catched!\n"); 128 } 129 130 } 131 /*}}}*/ 132 template class TPZRestoreClass<AdaptiveMeshRefinement,13829430>;/*{{{*/ 133 /*}}}*/ 134 void AdaptiveMeshRefinement::Write(TPZStream &buf,int withclassid){/*{{{*/ 135 136 try 137 { 138 /* Write context (this class) class ID*/ 139 TPZSaveable::Write(buf,withclassid); 140 /* Write this class id*/ 141 int classid = this->ClassId(); 142 buf.Write(&classid,1); 143 /* Write simple attributes */ 144 buf.Write(&this->levelmax,1); 145 buf.Write(&this->elementswidth,1); 146 buf.Write(&this->regionlevel1,1); 147 buf.Write(&this->regionlevelmax,1); 148 /* Write vector attributes*/ 149 int size=this->sid2index.size(); 150 int* psid2index=&this->sid2index[0]; 151 buf.Write(&size,1);//vector size 152 buf.Write(psid2index,this->sid2index.size()); 153 /* Write the geometric mesh*/ 154 this->fathermesh->Write(buf,this->ClassId()); 155 this->currentmesh->Write(buf,this->ClassId()); 156 } 157 catch(const std::exception& e) 158 { 159 _error_("AdaptiveMeshRefinement::Write: Exception catched!\n"); 160 } 80 this->index2sid.clear(); 81 this->specialelementsindex.clear(); 161 82 } 162 83 /*}}}*/ … … 176 97 if(numberofelements!=this->sid2index.size()) _error_("Impossible to execute refinement: sid2index.size is not equal to numberofelements!\n"); 177 98 178 /*Execute the refinement.*/ 179 this->RefinementProcess(amr_verbose,partiallyfloatedelements,masklevelset,deviatorictensorerror,thicknesserror); 180 181 /*Get new geometric mesh in ISSM data structure*/ 182 this->GetMesh(newnumberofvertices,newnumberofelements,x,y,elementslist); 183 184 /*Verify the new geometry*/ 185 this->CheckMesh(newnumberofvertices,newnumberofelements,this->elementswidth,x,y,elementslist); 186 187 } 188 /*}}}*/ 189 void AdaptiveMeshRefinement::RefinementProcess(bool &amr_verbose,double* partiallyfloatedelements,double* masklevelset, 190 double* deviatorictensorerror,double* thicknesserror){/*{{{*/ 191 192 if(amr_verbose) _printf_("\n\trefinement process started (level max = " << this->levelmax << ")\n"); 193 194 /*Intermediaries*/ 195 TPZGeoMesh* nohangingnodesmesh=NULL; 99 /*Combine the fields*/ 196 100 double mean_mask = 0; 197 101 double mean_tauerror = 0; 198 102 double mean_Herror = 0; 199 double group_error = 0; 200 103 int* fmask = xNew<int>(numberofelements); 104 int* ftauerror = xNew<int>(numberofelements); 105 int* fHerror = xNew<int>(numberofelements); 106 int* phi = xNew<int>(numberofelements); 201 107 /*Calculate mean values*/ 202 108 for(int i=0;i<this->sid2index.size();i++){ … … 208 114 mean_tauerror /= this->sid2index.size(); 209 115 mean_Herror /= this->sid2index.size(); 210 116 /*Filter to thickness*/ 117 for(int i=0;i<this->sid2index.size();i++){ 118 fmask[i]=0; 119 if(thicknesserror[i]>mean_Herror) fmask[i]=1; 120 } 121 /*Filter to tau*/ 122 for(int i=0;i<this->sid2index.size();i++){ 123 ftauerror[i]=0; 124 if(deviatorictensorerror[i]>mean_tauerror) ftauerror[i]=1; 125 } 126 /*Sum*/ 127 for(int i=0;i<this->sid2index.size();i++){ 128 phi[i]=ftauerror[i]+fmask[i]; 129 } 130 131 /*Execute the refinement.*/ 132 this->RefinementProcess(amr_verbose,partiallyfloatedelements,masklevelset,deviatorictensorerror,thicknesserror); 133 134 /*Get new geometric mesh in ISSM data structure*/ 135 this->GetMesh(newnumberofvertices,newnumberofelements,x,y,elementslist); 136 137 std::ofstream file("/home/santos/mesh.vtk"); 138 TPZVTKGeoMesh::PrintGMeshVTK(this->currentmesh,file); 139 140 /*Verify the new geometry*/ 141 this->CheckMesh(newnumberofvertices,newnumberofelements,this->elementswidth,x,y,elementslist); 142 143 } 144 /*}}}*/ 145 void AdaptiveMeshRefinement::RefinementProcess(bool &amr_verbose,double* partiallyfloatedelements,double* masklevelset,double* deviatorictensorerror,double* thicknesserror){/*{{{*/ 146 147 if(amr_verbose) _printf_("\n\trefinement process started (level max = " << this->levelmax << ")\n"); 148 149 /*Intermediaries*/ 150 TPZGeoMesh* nohangingnodesmesh=NULL; 151 152 //itapopo 153 this->step++; 154 155 double mean_mask = 0; 156 double mean_tauerror = 0; 157 double mean_Herror = 0; 158 double group_Herror = 0; 159 int index,sid; 160 std::vector<int> specialfatherindex; specialfatherindex.clear(); 161 162 /*Calculate mean values*/ 163 for(int i=0;i<this->sid2index.size();i++){ 164 mean_mask += masklevelset[i]; 165 mean_tauerror += deviatorictensorerror[i]; 166 mean_Herror += thicknesserror[i]; 167 } 168 mean_mask /= this->sid2index.size(); 169 mean_tauerror /= this->sid2index.size(); 170 mean_Herror /= this->sid2index.size(); 171 172 if(amr_verbose) _printf_("\t\tdeal with special elements...\n"); 173 /*Deal with special elements*/ 174 for(int i=0;i<this->specialelementsindex.size();i++){ 175 if(this->specialelementsindex[i]==-1) continue; 176 /*Get special element and verify*/ 177 TPZGeoEl* geoel=this->currentmesh->Element(this->specialelementsindex[i]); 178 if(!geoel)_error_("special element (sid) "<<i<<" is null!\n"); 179 if(geoel->HasSubElement())_error_("special element (sid) "<<i<<" has "<<geoel->NSubElements()<<" subelements!\n"); 180 if(geoel->MaterialId()!=this->GetElemMaterialID()) _error_("geoel->MaterialId is not GetElemMaterialID!\n"); 181 if(!geoel->Father())_error_("father of special element (sid) "<<i<<" is null!\n"); 182 183 /*Get element's siblings and verify*/ 184 TPZGeoEl* father=geoel->Father(); 185 TPZVec<TPZGeoEl *> siblings; 186 father->GetHigherSubElements(siblings); 187 std::vector<int> sidvec; sidvec.resize(siblings.size()); 188 for (int j=0;j<siblings.size();j++){ 189 if(!siblings[j]) _error_("special element (sid) "<<i<<" has a null siblings null!\n"); 190 sidvec[j]=this->index2sid[siblings[j]->Index()]; 191 } 192 193 /*Now, reset the data strucure and verify if the siblings should be deleted*/ 194 if(siblings.size()<4){ 195 /*Reset subelements in the father*/ 196 father->ResetSubElements(); 197 }else{ 198 if(siblings.size()!=4) _error_("element (index) "<<father->Index()<<" has "<<father->NSubElements()<<" subelements!\n"); 199 } 200 for (int j=0;j<siblings.size();j++){ 201 for(int k=0;k<this->specialelementsindex.size();k++){ 202 if(this->specialelementsindex[k]==siblings[j]->Index()){ 203 index = siblings[j]->Index(); 204 if(index<0) _error_("index is null!\n"); 205 sid = this->index2sid[index]; 206 if(sid<0) _error_("sid is null!\n"); 207 this->specialelementsindex[k] = -1; 208 this->index2sid[index] = -1; 209 this->sid2index[sid] = -1; 210 } 211 } 212 if(siblings.size()<4){ 213 /*Ok, the special element can be deleted*/ 214 siblings[j]->ResetSubElements(); 215 this->currentmesh->DeleteElement(siblings[j],siblings[j]->Index()); 216 } 217 } 218 219 /*Now, verify if the father should be refined with uniform pattern (smoother)*/ 220 this->smooth_frequency=10000000; 221 if(siblings.size()==3 || this->step%this->smooth_frequency==0){//it keeps the mesh with uniform elements 222 /*Father has uniform subelements now*/ 223 TPZVec<TPZGeoEl *> sons; 224 father->Divide(sons); 225 this->smooth_frequency=this->step; 226 }else{ 227 specialfatherindex.push_back(father->Index()); 228 } 229 if(this->specialelementsindex[i]!=-1) _error_("special element "<<i<<" was not deleted!\n"); 230 } 231 this->currentmesh->BuildConnectivity(); 232 211 233 if(amr_verbose) _printf_("\t\tuniform refinement...\n"); 234 /*Deal with uniform elemnts*/ 212 235 for(int i=0;i<this->sid2index.size();i++){ 236 if(this->sid2index[i]==-1) continue; 237 /*Get element and verify*/ 213 238 TPZGeoEl* geoel=this->currentmesh->Element(this->sid2index[i]); 214 if(geoel->HasSubElement()) _error_("Impossible to refine: geoel (index) "<<this->sid2index[i]<<" has subelements!\n"); 215 if(geoel->MaterialId()!=this->GetElemMaterialID()) _error_("Impossible to refine: geoel->MaterialId is not GetElemMaterialID!\n"); 216 /*Refine*/ 217 if(thicknesserror[i]>mean_Herror){ 239 if(geoel->HasSubElement()) _error_("element (sid) "<<i<<" has "<<geoel->NSubElements()<<" subelements!\n"); 240 if(geoel->MaterialId()!=this->GetElemMaterialID()) _error_("geoel->MaterialId is not GetElemMaterialID!\n"); 241 242 /*Refine process*/ 243 if(thicknesserror[i]>mean_Herror) 244 { 245 int count=0; 246 TPZGeoEl* father=geoel->Father(); 247 if(father){ 248 for(int j=3;j<6;j++){ 249 index=father->Neighbour(j).Element()->Index(); 250 for(int k=0;k<specialfatherindex.size();k++) if(specialfatherindex[k]==index) count++; 251 } 252 } 218 253 TPZVec<TPZGeoEl *> sons; 219 if(geoel->Level()<this->levelmax ) geoel->Divide(sons);254 if(geoel->Level()<this->levelmax && count==0) geoel->Divide(sons); 220 255 } 221 else if(geoel->Level()>0){ /*try to unrefine*/ 222 TPZVec<TPZGeoEl *> sons; 223 geoel->Father()->GetHigherSubElements(sons); 224 group_error=0; 225 for(int j=0;j<sons.size();j++){ 226 sons[j]->Index(); 256 else if(geoel->Level()>0) 257 {/*Unrefine process*/ 258 259 /*Get siblings and verify*/ 260 TPZVec<TPZGeoEl *> siblings; 261 geoel->Father()->GetHigherSubElements(siblings); 262 //if(siblings.size()<4) _error_("Impossible to refine: geoel (index) "<<this->sid2index[i]<<" has less than 3 siblings!\n"); 263 if(siblings.size()>4) continue;//Father has more then 4 sons, this group should not be unrefined. 264 265 /*Compute the error of the group*/ 266 group_Herror=0; 267 for(int j=0;j<siblings.size();j++){ 268 index = siblings[j]->Index(); 269 sid = this->index2sid[index]; 270 if(sid==-1) continue; 271 group_Herror+=thicknesserror[sid]; 227 272 } 228 } 229 } 273 /*Verify if this group should be unrefined*/ 274 if(group_Herror>0 && group_Herror<0*mean_Herror){ //itapopo 275 /*Reset subelements in the father*/ 276 this->currentmesh->Element(geoel->Father()->Index())->ResetSubElements(); 277 /*Delete the elements and set their indexes in the index2sid and sid2index*/ 278 for (int j=0;j<siblings.size();j++){ 279 index = siblings[j]->Index(); 280 sid = this->index2sid[index]; 281 this->index2sid[index]=-1; 282 if(sid!=-1) this->sid2index[sid]=-1; 283 this->currentmesh->DeleteElement(siblings[j],siblings[j]->Index()); 284 }//for j 285 }//if 286 }/*Unrefine process*/ 287 }//for i 230 288 this->currentmesh->BuildConnectivity(); 231 289 232 290 if(amr_verbose) _printf_("\t\trefine to avoid hanging nodes...\n"); 233 291 this->RefineMeshToAvoidHangingNodes(this->currentmesh); 234 this->currentmesh->BuildConnectivity();235 292 236 293 //nohangingnodesmesh = this->CreateRefPatternMesh(newmesh); itapopo tentar otimizar 237 294 238 295 if(amr_verbose) _printf_("\trefinement process done!\n"); 239 296 } … … 257 314 void AdaptiveMeshRefinement::RefineMeshToAvoidHangingNodes(TPZGeoMesh *gmesh){/*{{{*/ 258 315 259 /*Refine elements to avoid hanging nodes: non-uniform refinement*/ 316 /*Now, insert special elements to avoid hanging nodes*/ 317 this->specialelementsindex.clear(); 260 318 const int NElem = gmesh->NElements(); 261 319 for(int i=0;i<NElem;i++){ … … 268 326 TPZAutoPointer<TPZRefPattern> refp=TPZRefPatternTools::PerfectMatchRefPattern(geoel); 269 327 if(refp){ 328 /*Non-uniform refinement*/ 270 329 TPZVec<TPZGeoEl *> Sons; 271 330 geoel->SetRefPattern(refp); 272 331 geoel->Divide(Sons); 273 } 274 } 275 gmesh->BuildConnectivity(); 276 332 /*Keep the index of the special elements*/ 333 for(int j=0;j<Sons.size();j++) this->specialelementsindex.push_back(Sons[j]->Index()); 334 } 335 } 336 gmesh->BuildConnectivity(); 277 337 } 278 338 /*}}}*/ … … 281 341 /* IMPORTANT! pelements are in Matlab indexing 282 342 NEOPZ works only in C indexing. 283 This method cleans up and updated the this->sid2index 284 and fills in it with the new mesh.343 This method cleans up and updated the this->sid2index 344 and this->index2sid and fills in it with the new mesh. 285 345 Avoid to call this method before Refinement Process.*/ 286 346 … … 295 355 TPZGeoEl* geoel = NULL; 296 356 long* vertex_index2sid = xNew<long>(gmesh->NNodes()); 357 this->index2sid.clear(); this->index2sid.resize(gmesh->NElements()); 297 358 this->sid2index.clear(); 298 359 … … 307 368 /*Fill in the vertex_index2sid vector with non usual index value*/ 308 369 for(int i=0;i<gmesh->NNodes();i++) vertex_index2sid[i]=-1; 370 371 /*Fill in the this->index2sid vector with non usual index value*/ 372 for(int i=0;i<gmesh->NElements();i++) this->index2sid[i]=-1; 309 373 310 374 /*Get elements without sons and fill in the vertex_index2sid with used vertices (indexes) */ … … 316 380 if(geoel->MaterialId() != this->GetElemMaterialID()) continue; 317 381 this->sid2index.push_back(i);//keep the element index 382 this->index2sid[i]=this->sid2index.size()-1;//keep the element sid 318 383 for(int j=0;j<this->elementswidth;j++){ 319 384 nodeindex=geoel->NodeIndex(j); … … 347 412 newelements[i*this->elementswidth+j]=(int)sid+1;//C to Matlab indexing 348 413 } 349 } 350 414 /*Verify the Jacobian determinant. If detJ<0, swap the 2 first postions: 415 a -> b 416 b -> a */ 417 double detJ,xa,xb,xc,ya,yb,yc; 418 int a,b,c; 419 420 a=newelements[i*this->elementswidth+0]-1; 421 b=newelements[i*this->elementswidth+1]-1; 422 c=newelements[i*this->elementswidth+2]-1; 423 424 xa=newmeshX[a]; ya=newmeshY[a]; 425 xb=newmeshX[b]; yb=newmeshY[b]; 426 xc=newmeshX[c]; yc=newmeshY[c]; 427 428 detJ=(xb-xa)*(yc-ya)-(xc-xa)*(yb-ya); 429 430 /*verify and swap, if necessary*/ 431 if(detJ<0) { 432 newelements[i*this->elementswidth+0]=b+1;//a->b 433 newelements[i*this->elementswidth+1]=a+1;//b->a 434 } 435 } 436 351 437 /*Setting outputs*/ 352 438 nvertices = nconformvertices; … … 358 444 /*Cleanup*/ 359 445 xDelete<long>(vertex_index2sid); 360 361 446 } 362 447 /*}}}*/ -
issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h
r21806 r21919 12 12 /*REAL and STATE definitions, NeoPZ variables itapopo should be read by NeoPZ's config.h*/ 13 13 #ifndef REFPATTERNDIR 14 # 14 #define REFPATTERNDIR "/home/santos/trunk-jpl/externalpackages/neopz/install/include/refpatterns" 15 15 #endif 16 16 … … 24 24 25 25 #include <pzreal.h> 26 #include <pzsave.h>27 26 #include <pzgmesh.h> 28 27 #include <pzvec.h> … … 36 35 #include <TPZGeoElement.h> 37 36 #include <pzreftriangle.h> 37 #include <pzgeotriangle.h> 38 38 #include <tpzgeoelrefpattern.h> 39 #include <TPZVTKGeoMesh.h> 40 41 #include "../shared/shared.h" 42 39 43 /*}}}*/ 40 44 41 class AdaptiveMeshRefinement : public TPZSaveable{45 class AdaptiveMeshRefinement{ 42 46 43 47 public: … … 49 53 AdaptiveMeshRefinement & operator= (const AdaptiveMeshRefinement &cp); 50 54 virtual ~AdaptiveMeshRefinement(); 51 52 /*Savable methods*/53 virtual int ClassId() const;54 virtual void Read(TPZStream &buf,void *context);55 virtual void Write(TPZStream &buf,int withclassid);56 55 57 56 /*General methods*/ … … 70 69 private: 71 70 /*Private attributes*/ 72 71 int elementswidth; // Geometric nodes for element: 3 == Tria, 4 == Tetra, 6 == Penta 73 72 int levelmax; // Max level of refinement 74 73 double regionlevel1; // Region which will be refined with level 1 75 74 double regionlevelmax; // Region which will be refined with level max 75 int step; //itapopo testando 76 int smooth_frequency; //itapopo testando 76 77 std::vector<int> sid2index; // Vector that keeps index of PZGeoMesh elements used in the ISSM mesh (sid) 78 std::vector<int> index2sid; // Vector that keeps sid of issm mesh elements used in the neopz mesh (index) 79 std::vector<int> specialelementsindex; // Vector that keeps index of the special elements (created to avoid haning nodes) 77 80 TPZGeoMesh *fathermesh; // Father Mesh is the entire mesh without refinement 78 81 TPZGeoMesh *currentmesh; // Current Mesh is the refined mesh -
issm/trunk-jpl/src/c/classes/AmrBamg.cpp
r21864 r21919 17 17 18 18 /*Constructor, copy, clean up and destructor*/ 19 AmrBamg::AmrBamg(IssmDouble hmin, IssmDouble hmax,int fieldenum_in,IssmDouble err){/*{{{*/ 19 AmrBamg::AmrBamg(IssmDouble hmin,IssmDouble hmax,int fieldenum_in,IssmDouble err_in,int keepmetric_in, 20 IssmDouble groundingline_resolution_in,IssmDouble groundingline_distance_in, 21 IssmDouble icefront_resolution_in,IssmDouble icefront_distance_in, 22 IssmDouble thicknesserror_resolution_in,IssmDouble thicknesserror_threshold_in, 23 IssmDouble deviatoricerror_resolution_in,IssmDouble deviatoricerror_threshold_in){/*{{{*/ 20 24 21 this->fieldenum = fieldenum_in; 22 this->geometry = NULL; 23 this->fathermesh = NULL; 24 this->previousmesh = NULL; 25 this->fieldenum = fieldenum_in; 26 this->keepmetric = keepmetric_in; 27 this->groundingline_resolution = groundingline_resolution_in; 28 this->groundingline_distance = groundingline_distance_in; 29 this->icefront_resolution = icefront_resolution_in; 30 this->icefront_distance = icefront_distance_in; 31 this->thicknesserror_resolution = thicknesserror_resolution_in; 32 this->thicknesserror_threshold = thicknesserror_threshold_in; 33 this->deviatoricerror_resolution = deviatoricerror_resolution_in; 34 this->deviatoricerror_threshold = deviatoricerror_threshold_in; 35 this->geometry = NULL; 36 this->fathermesh = NULL; 37 this->previousmesh = NULL; 25 38 26 39 /*Only initialize options for now (same as bamg.m)*/ … … 49 62 this->options->errSize[0]=1; 50 63 this->options->errSize[1]=1; 51 this->options->err[0] = err ;64 this->options->err[0] = err_in; 52 65 } 53 66 /*}}}*/ … … 83 96 delete Th; 84 97 }/*}}}*/ 85 void AmrBamg::ExecuteRefinementBamg(IssmDouble* field, int* pnewnumberofvertices,int *pnewnumberofelements,IssmDouble** px,IssmDouble** py,IssmDouble** pz,int** pelementslist){/*{{{*/98 void AmrBamg::ExecuteRefinementBamg(IssmDouble* field,IssmDouble* hmaxVertices,int* pnewnumberofvertices,int *pnewnumberofelements,IssmDouble** px,IssmDouble** py,IssmDouble** pz,int** pelementslist){/*{{{*/ 86 99 87 100 /*Intermediaries*/ … … 93 106 _assert_(this->options); 94 107 _assert_(this->fathermesh); 95 _assert_(field );108 _assert_(field || hmaxVertices);//at least one is necessary 96 109 97 110 /*Prepare field for metric*/ 98 this->options->field = field; 111 this->options->field = field; 112 this->options->hmaxVertices = hmaxVertices; 99 113 100 114 /*remesh*/ 101 115 if(this->previousmesh){ 102 this->options->fieldSize[0] = this->previousmesh->VerticesSize[0]; 103 this->options->fieldSize[1] = 1; 116 this->options->fieldSize[0] = this->previousmesh->VerticesSize[0]; 117 this->options->fieldSize[1] = 1; 118 this->options->hmaxVerticesSize[0] = this->previousmesh->VerticesSize[0]; 119 this->options->hmaxVerticesSize[1] = 1; 104 120 Bamgx(meshout,geomout,this->previousmesh,this->geometry,this->options); 105 121 } 106 122 else{ 107 this->options->fieldSize[0] = this->fathermesh->VerticesSize[0]; 108 this->options->fieldSize[1] = 1; 123 this->options->fieldSize[0] = this->fathermesh->VerticesSize[0]; 124 this->options->fieldSize[1] = 1; 125 this->options->hmaxVerticesSize[0] = this->fathermesh->VerticesSize[0]; 126 this->options->hmaxVerticesSize[1] = 1; 109 127 Bamgx(meshout,geomout,this->fathermesh,this->geometry,this->options); 110 128 } 111 129 112 /*remove field for memory management (FemModel is taking care of deleting it)*/130 /*remove field and hmaxVertices for memory management (FemModel is taking care of deleting it)*/ 113 131 this->options->field = NULL; 114 132 this->options->fieldSize[0] = 0; 115 133 this->options->fieldSize[1] = 0; 134 this->options->hmaxVertices = NULL; 135 this->options->hmaxVerticesSize[0] = 0; 136 this->options->hmaxVerticesSize[1] = 0; 137 138 /*verify if the metric will be reseted or not*/ 139 if(!this->keepmetric){ 140 if(this->options->metric) xDelete<IssmDouble>(this->options->metric); 141 this->options->metricSize[0] = 0; 142 this->options->metricSize[1] = 0; 143 } 116 144 117 145 /*Change previous mesh*/ -
issm/trunk-jpl/src/c/classes/AmrBamg.h
r21812 r21919 14 14 public: 15 15 int fieldenum; 16 int keepmetric; 17 IssmDouble groundingline_resolution; 18 IssmDouble groundingline_distance; 19 IssmDouble icefront_resolution; 20 IssmDouble icefront_distance; 21 IssmDouble thicknesserror_resolution; 22 IssmDouble thicknesserror_threshold; 23 IssmDouble deviatoricerror_resolution; 24 IssmDouble deviatoricerror_threshold; 16 25 17 26 /* Constructor, destructor etc*/ 18 AmrBamg(IssmDouble hmin, IssmDouble hmax,int fieldenum_in,IssmDouble err_in); 27 AmrBamg(IssmDouble hmin,IssmDouble hmax,int fieldenum_in,IssmDouble err_in,int keepmetric_in, 28 IssmDouble groundingline_resolution_in,IssmDouble groundingline_distance_in, 29 IssmDouble icefront_resolution_in,IssmDouble icefront_distance_in, 30 IssmDouble thicknesserror_resolution_in,IssmDouble thicknesserror_threshold_in, 31 IssmDouble deviatoricerror_resolution_in,IssmDouble deviatoricerror_threshold_in); 19 32 ~AmrBamg(); 20 33 21 34 /*General methods*/ 22 35 void Initialize(int* elements,IssmDouble* x,IssmDouble* y,int numberofvertices,int numberofelements); 23 void ExecuteRefinementBamg(IssmDouble* field, int* pnewnumberofvertices,int *pnewnumberofelements,IssmDouble** px,IssmDouble** py,IssmDouble** pz,int** pelementslist);36 void ExecuteRefinementBamg(IssmDouble* field,IssmDouble* hmaxVertices,int* pnewnumberofvertices,int *pnewnumberofelements,IssmDouble** px,IssmDouble** py,IssmDouble** pz,int** pelementslist); 24 37 25 38 private: -
issm/trunk-jpl/src/c/classes/Elements/Element.cpp
r21917 r21919 1711 1711 name==LevelsetfunctionSlopeYEnum || 1712 1712 name==LevelsetfunctionPicardEnum || 1713 //name==CalvingCalvingrateEnum || 1713 name==CalvingCalvingrateEnum || 1714 name==CalvingMeltingrateEnum || 1714 1715 name==GradientEnum || 1715 1716 name==OldGradientEnum || -
issm/trunk-jpl/src/c/classes/FemModel.cpp
r21915 r21919 2303 2303 #endif 2304 2304 2305 2305 #if defined(_HAVE_BAMG_) && !defined(_HAVE_ADOLC_) 2306 2306 case AmrBamgEnum: this->ReMeshBamg(&newnumberofvertices,&newnumberofelements,&newx,&newy,&newz,&newelementslist); break; 2307 2307 #endif 2308 2308 2309 2309 default: _error_("not implemented yet"); … … 2720 2720 void FemModel::WriteMeshInResults(void){/*{{{*/ 2721 2721 2722 //itapopo2723 2722 #if defined(_HAVE_NEOPZ_) && !defined(_HAVE_ADOLC_) 2724 2723 this->WriteErrorEstimatorsInResults(); 2725 2724 #endif 2726 //itapopo2727 2725 2728 2726 int step = -1; … … 2754 2752 this->results->AddResult(new GenericExternalResult<IssmDouble*>(this->results->Size()+1,MeshYEnum, 2755 2753 y,numberofvertices,1,step,time)); 2756 2757 //itapopo2758 #if defined(_HAVE_NEOPZ_) && !defined(_HAVE_ADOLC_)2759 if(IssmComm::GetRank()==0){2760 TPZFileStream fstr;2761 std::stringstream ss;2762 int frictionlaw,levelmax;2763 this->parameters->FindParam(&frictionlaw,FrictionLawEnum);2764 this->parameters->FindParam(&levelmax,AmrLevelMaxEnum);2765 ss << levelmax;2766 if(frictionlaw==1){2767 ss << "_viscous/amr.txt";2768 }else if(frictionlaw==7){2769 ss << "_tsai/amr.txt";2770 }else{2771 _error_("friction law not supported here.");2772 }2773 std::string AMRfile = "/home/santos/Misomip2/L" + ss.str();2774 fstr.OpenWrite(AMRfile.c_str());2775 int withclassid = 1;2776 this->amr->Write(fstr,withclassid);2777 }2778 #endif2779 //itapopo2780 2754 2781 2755 /*Cleanup*/ … … 3614 3588 } 3615 3589 /*}}}*/ 3590 void FemModel::GetElementCenterCoordinates(IssmDouble** pxc,IssmDouble** pyc){/*{{{*/ 3591 3592 /*Intermediaries*/ 3593 int elementswidth = this->GetElementsWidth(); 3594 int numberofelements = this->elements->NumberOfElements(); 3595 int* elem_vertices = xNew<int>(elementswidth); 3596 Vector<IssmDouble>* vxc = new Vector<IssmDouble>(numberofelements); 3597 Vector<IssmDouble>* vyc = new Vector<IssmDouble>(numberofelements); 3598 IssmDouble *x = NULL; 3599 IssmDouble *y = NULL; 3600 IssmDouble *z = NULL; 3601 3602 /*Get vertices coordinates*/ 3603 VertexCoordinatesx(&x,&y,&z,this->vertices,false) ; 3604 3605 /*Insert the element center coordinates*/ 3606 for(int i=0;i<this->elements->Size();i++){ 3607 Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i)); 3608 element->GetVerticesSidList(elem_vertices); 3609 int sid = element->Sid(); 3610 vxc->SetValue(sid,(x[elem_vertices[0]]+x[elem_vertices[1]]+x[elem_vertices[2]])/3.,INS_VAL); 3611 vyc->SetValue(sid,(y[elem_vertices[0]]+y[elem_vertices[1]]+y[elem_vertices[2]])/3.,INS_VAL); 3612 } 3613 3614 /*Assemble*/ 3615 vxc->Assemble(); 3616 vyc->Assemble(); 3617 3618 /*Serialize and set output*/ 3619 (*pxc)=vxc->ToMPISerial(); 3620 (*pyc)=vyc->ToMPISerial(); 3621 3622 /*Cleanup*/ 3623 xDelete<IssmDouble>(x); 3624 xDelete<IssmDouble>(y); 3625 xDelete<IssmDouble>(z); 3626 xDelete<int>(elem_vertices); 3627 delete vxc; 3628 delete vyc; 3629 } 3630 /*}}}*/ 3616 3631 #endif 3617 3632 … … 4336 4351 int my_rank = IssmComm::GetRank(); 4337 4352 4353 /*Intermediaries*/ 4354 int numberofvertices = this->vertices->NumberOfVertices(); 4355 IssmDouble* vector_serial = NULL; 4356 IssmDouble* hmaxvertices_serial = NULL; 4357 Vector<IssmDouble> *vector = NULL; 4358 4338 4359 /*Get vector to create metric*/ 4339 int numberofvertices = this->vertices->NumberOfVertices(); 4340 Vector<IssmDouble> *vector = NULL; 4341 GetVectorFromInputsx(&vector,this,this->amrbamg->fieldenum,VertexSIdEnum); 4342 vector->Assemble(); 4343 IssmDouble* vector_serial = vector->ToMPISerial(); 4360 if(this->amrbamg->fieldenum!=NoneEnum){ 4361 GetVectorFromInputsx(&vector,this,this->amrbamg->fieldenum,VertexSIdEnum); 4362 vector->Assemble(); 4363 vector_serial = vector->ToMPISerial(); 4364 } 4365 4366 /*Get hmaxVertices to create metric*/ 4367 if(this->amrbamg->groundingline_distance>0||this->amrbamg->icefront_distance>0|| 4368 this->amrbamg->thicknesserror_threshold>0||this->amrbamg->deviatoricerror_threshold>0){ 4369 /*Initialize hmaxvertices with NAN*/ 4370 hmaxvertices_serial=xNew<IssmDouble>(numberofvertices); 4371 for(int i=0;i<numberofvertices;i++) hmaxvertices_serial[i]=NAN; 4372 /*Fill hmaxvertices*/ 4373 if(this->amrbamg->groundingline_distance>0) this->GethmaxVerticesFromZeroLevelSetDistance(hmaxvertices_serial,MaskGroundediceLevelsetEnum); 4374 if(this->amrbamg->icefront_distance>0) this->GethmaxVerticesFromZeroLevelSetDistance(hmaxvertices_serial,MaskIceLevelsetEnum); 4375 if(this->amrbamg->thicknesserror_threshold>0) this->GethmaxVerticesFromEstimators(hmaxvertices_serial,ThicknessErrorEstimatorEnum); 4376 if(this->amrbamg->deviatoricerror_threshold>0) this->GethmaxVerticesFromEstimators(hmaxvertices_serial,DeviatoricStressErrorEstimatorEnum); 4377 } 4378 4379 if(my_rank==0){ 4380 this->amrbamg->ExecuteRefinementBamg(vector_serial,hmaxvertices_serial,&newnumberofvertices,&newnumberofelements,&newx,&newy,&newz,&newelementslist); 4381 if(newnumberofvertices<=0 || newnumberofelements<=0) _error_("Error in the refinement process."); 4382 } 4383 4384 /*Cleanup*/ 4385 xDelete<IssmDouble>(vector_serial); 4386 xDelete<IssmDouble>(hmaxvertices_serial); 4344 4387 delete vector; 4345 4346 if(my_rank==0){4347 this->amrbamg->ExecuteRefinementBamg(vector_serial,&newnumberofvertices,&newnumberofelements,&newx,&newy,&newz,&newelementslist);4348 if(newnumberofvertices<=0 || newnumberofelements<=0) _error_("Error in the refinement process.");4349 }4350 4351 xDelete<IssmDouble>(vector_serial);4352 4388 4353 4389 /*Send new mesh to others CPU*/ … … 4379 4415 int numberofvertices = this->vertices->NumberOfVertices(); 4380 4416 int numberofelements = this->elements->NumberOfElements(); 4381 int numberofsegments = 0; //used on matlab4382 4417 IssmDouble* x = NULL; 4383 4418 IssmDouble* y = NULL; 4384 4419 IssmDouble* z = NULL; 4385 4420 int* elements = NULL; 4386 int elementswidth = this->GetElementsWidth(); //just tria elements in this version. Itapopo:4387 4421 IssmDouble hmin,hmax,err; 4388 int fieldenum; 4422 IssmDouble groundingline_resolution,groundingline_distance,icefront_resolution,icefront_distance; 4423 IssmDouble thicknesserror_resolution,thicknesserror_threshold,deviatoricerror_resolution,deviatoricerror_threshold; 4424 int fieldenum,keepmetric; 4389 4425 4390 4426 /*Get rank*/ … … 4402 4438 this->parameters->FindParam(&fieldenum,AmrFieldEnum); 4403 4439 this->parameters->FindParam(&err,AmrErrEnum); 4440 this->parameters->FindParam(&keepmetric,AmrKeepMetricEnum); 4441 this->parameters->FindParam(&groundingline_resolution,AmrGroundingLineResolutionEnum); 4442 this->parameters->FindParam(&groundingline_distance,AmrGroundingLineDistanceEnum); 4443 this->parameters->FindParam(&icefront_resolution,AmrIceFrontResolutionEnum); 4444 this->parameters->FindParam(&icefront_distance,AmrIceFrontDistanceEnum); 4445 this->parameters->FindParam(&thicknesserror_resolution,AmrThicknessErrorResolutionEnum); 4446 this->parameters->FindParam(&thicknesserror_threshold,AmrThicknessErrorThresholdEnum); 4447 this->parameters->FindParam(&deviatoricerror_resolution,AmrDeviatoricErrorResolutionEnum); 4448 this->parameters->FindParam(&deviatoricerror_threshold,AmrDeviatoricErrorThresholdEnum); 4404 4449 4405 4450 /*Create bamg data structures for bamg*/ 4406 this->amrbamg = new AmrBamg(hmin,hmax,fieldenum,err); 4451 this->amrbamg = new AmrBamg(hmin,hmax,fieldenum,err,keepmetric, 4452 groundingline_resolution,groundingline_distance, 4453 icefront_resolution,icefront_distance, 4454 thicknesserror_resolution,thicknesserror_threshold, 4455 deviatoricerror_resolution,deviatoricerror_threshold); 4407 4456 4408 4457 /*Re-create original mesh and put it in bamg structure (only cpu 0)*/ … … 4416 4465 xDelete<IssmDouble>(z); 4417 4466 xDelete<int>(elements); 4467 } 4468 /*}}}*/ 4469 void FemModel::GethmaxVerticesFromZeroLevelSetDistance(IssmDouble* hmaxvertices,int levelset_type){/*{{{*/ 4470 4471 if(!hmaxvertices) _error_("hmaxvertices is NULL!\n"); 4472 4473 /*Intermediaries*/ 4474 int numberofvertices = this->vertices->NumberOfVertices(); 4475 IssmDouble* verticedistance = NULL; 4476 IssmDouble threshold,resolution; 4477 4478 switch(levelset_type){ 4479 case MaskGroundediceLevelsetEnum: 4480 threshold = this->amrbamg->groundingline_distance; 4481 resolution = this->amrbamg->groundingline_resolution; 4482 break; 4483 case MaskIceLevelsetEnum: 4484 threshold = this->amrbamg->icefront_distance; 4485 resolution = this->amrbamg->icefront_resolution; 4486 break; 4487 default: _error_("not implemented yet"); 4488 } 4489 4490 /*Get vertice distance to zero level set points*/ 4491 this->GetVerticeDistanceToZeroLevelSet(&verticedistance,levelset_type); 4492 if(!verticedistance) _error_("verticedistance is NULL!\n"); 4493 4494 /*Fill hmaxVertices*/ 4495 for(int i=0;i<numberofvertices;i++){ 4496 if(verticedistance[i]<threshold){ 4497 if(xIsNan<IssmDouble>(hmaxvertices[i])) hmaxvertices[i]=resolution; 4498 else hmaxvertices[i]=min(resolution,hmaxvertices[i]); 4499 } 4500 } 4501 4502 /*Cleanup*/ 4503 xDelete<IssmDouble>(verticedistance); 4504 } 4505 /*}}}*/ 4506 void FemModel::GethmaxVerticesFromEstimators(IssmDouble* hmaxvertices,int errorestimator_type){/*{{{*/ 4507 4508 if(!hmaxvertices) _error_("hmaxvertices is NULL!\n"); 4509 4510 /*Intermediaries*/ 4511 int numberofelements = this->elements->NumberOfElements(); 4512 int numberofvertices = this->vertices->NumberOfVertices(); 4513 IssmDouble* error_elements = NULL; 4514 IssmDouble* error_vertices = NULL; 4515 IssmDouble *x = NULL; 4516 IssmDouble *y = NULL; 4517 IssmDouble *z = NULL; 4518 int *index = NULL; 4519 IssmDouble maxerror,threshold,resolution; 4520 4521 switch(errorestimator_type){ 4522 case ThicknessErrorEstimatorEnum: 4523 threshold = this->amrbamg->thicknesserror_threshold; 4524 resolution = this->amrbamg->thicknesserror_resolution; 4525 this->ThicknessZZErrorEstimator(&error_elements); 4526 break; 4527 case DeviatoricStressErrorEstimatorEnum: 4528 threshold = this->amrbamg->deviatoricerror_threshold; 4529 resolution = this->amrbamg->deviatoricerror_resolution; 4530 this->ZZErrorEstimator(&error_elements); 4531 break; 4532 default: _error_("not implemented yet"); 4533 } 4534 4535 if(!error_elements) _error_("error_elements is NULL!\n"); 4536 4537 /*Get mesh*/ 4538 this->GetMesh(this->vertices,this->elements,&x,&y,&z,&index); 4539 4540 /*Sum the estimators over the vertices*/ 4541 error_vertices=xNewZeroInit<IssmDouble>(numberofvertices); 4542 for(int i=0;i<numberofelements;i++){ 4543 for(int j=0;j<this->GetElementsWidth();j++){ 4544 int vid=index[i*this->GetElementsWidth()+j]-1;//Matlab to C indexing 4545 error_vertices[vid]+=error_elements[i]; 4546 } 4547 } 4548 4549 /*Find the max of the estimators (use error_elements)*/ 4550 maxerror=error_elements[0]; 4551 for(int i=0;i<numberofelements;i++) maxerror=max(maxerror,error_elements[i]); 4552 4553 /*Fill hmaxvertices*/ 4554 for(int i=0;i<numberofvertices;i++){ 4555 if(error_vertices[i]>threshold*maxerror){ 4556 if(xIsNan<IssmDouble>(hmaxvertices[i])) hmaxvertices[i]=resolution; 4557 else hmaxvertices[i]=min(resolution,hmaxvertices[i]); 4558 } 4559 } 4560 4561 /*Cleanup*/ 4562 xDelete<IssmDouble>(x); 4563 xDelete<IssmDouble>(y); 4564 xDelete<IssmDouble>(z); 4565 xDelete<int>(index); 4566 xDelete<IssmDouble>(error_elements); 4567 xDelete<IssmDouble>(error_vertices); 4568 } 4569 /*}}}*/ 4570 void FemModel::GetVerticeDistanceToZeroLevelSet(IssmDouble** pverticedistance,int levelset_type){/*{{{*/ 4571 4572 /*Here, "zero level set" means grounding line or ice front, depending on the level set type*/ 4573 /*pverticedistance is the minimal vertice distance to the grounding line or ice front*/ 4574 if(levelset_type!=MaskGroundediceLevelsetEnum && levelset_type!=MaskIceLevelsetEnum){ 4575 _error_("level set type not implemented yet!"); 4576 } 4577 4578 /*Output*/ 4579 IssmDouble* verticedistance; 4580 4581 /*Intermediaries*/ 4582 int elementswidth = this->GetElementsWidth(); 4583 int numberofvertices = this->vertices->NumberOfVertices(); 4584 int numberofelements = this->elements->NumberOfElements(); 4585 int* elem_vertices = xNew<int>(elementswidth); 4586 IssmDouble* levelset = xNew<IssmDouble>(elementswidth); 4587 IssmDouble* xp = NULL; 4588 IssmDouble* yp = NULL; 4589 IssmDouble* x = NULL; 4590 IssmDouble* y = NULL; 4591 IssmDouble* z = NULL; 4592 Vector<IssmDouble>* vx_zerolevelset = new Vector<IssmDouble>(numberofelements); 4593 Vector<IssmDouble>* vy_zerolevelset = new Vector<IssmDouble>(numberofelements); 4594 IssmDouble* x_zerolevelset = NULL; 4595 IssmDouble* y_zerolevelset = NULL; 4596 int count,sid,npoints; 4597 IssmDouble xcenter,ycenter,distance; 4598 4599 /*Get vertices coordinates*/ 4600 VertexCoordinatesx(&x,&y,&z,this->vertices,false) ; 4601 4602 /*Use the element center coordinate if level set is zero (grounding line or ice front), otherwise set NAN*/ 4603 for(int i=0;i<this->elements->Size();i++){ 4604 Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i)); 4605 element->GetInputListOnVertices(levelset,levelset_type); 4606 element->GetVerticesSidList(elem_vertices); 4607 sid = element->Sid(); 4608 xcenter = NAN; 4609 ycenter = NAN; 4610 Tria* tria = xDynamicCast<Tria*>(element); 4611 if(tria->IsIceInElement()){/*verify if there is ice in the element*/ 4612 if(levelset[0]*levelset[1]<0. || levelset[0]*levelset[2]<0. || 4613 abs(levelset[0]*levelset[1])<DBL_EPSILON || abs(levelset[0]*levelset[2])<DBL_EPSILON) { 4614 xcenter=(x[elem_vertices[0]]+x[elem_vertices[1]]+x[elem_vertices[2]])/3.; 4615 ycenter=(y[elem_vertices[0]]+y[elem_vertices[1]]+y[elem_vertices[2]])/3.; 4616 } 4617 } 4618 vx_zerolevelset->SetValue(sid,xcenter,INS_VAL); 4619 vy_zerolevelset->SetValue(sid,ycenter,INS_VAL); 4620 } 4621 /*Assemble and serialize*/ 4622 vx_zerolevelset->Assemble(); 4623 vy_zerolevelset->Assemble(); 4624 x_zerolevelset=vx_zerolevelset->ToMPISerial(); 4625 y_zerolevelset=vy_zerolevelset->ToMPISerial(); 4626 4627 /*keep just the element center coordinates with zero level set (compact the structure to save time in the next step)*/ 4628 count = 0; 4629 xp = xNewZeroInit<IssmDouble>(numberofelements); 4630 yp = xNewZeroInit<IssmDouble>(numberofelements); 4631 for(int i=0;i<numberofelements;i++){ 4632 if(!xIsNan<IssmDouble>(x_zerolevelset[i])){ 4633 xp[count]=x_zerolevelset[i]; 4634 yp[count]=y_zerolevelset[i]; 4635 count++; 4636 } 4637 } 4638 npoints=count; 4639 4640 /*Find the minimal vertice distance to the zero levelset (grounding line or ice front)*/ 4641 verticedistance=xNew<IssmDouble>(numberofvertices); 4642 for(int i=0;i<numberofvertices;i++){ 4643 verticedistance[i]=INFINITY; 4644 for(int j=0;j<npoints;j++){ 4645 distance=sqrt((x[i]-xp[j])*(x[i]-xp[j])+(y[i]-yp[j])*(y[i]-yp[j])); 4646 verticedistance[i]=min(distance,verticedistance[i]); 4647 } 4648 } 4649 4650 /*Assign the pointer*/ 4651 (*pverticedistance)=verticedistance; 4652 4653 /*Cleanup*/ 4654 xDelete<int>(elem_vertices); 4655 xDelete<IssmDouble>(levelset); 4656 xDelete<IssmDouble>(x_zerolevelset); 4657 xDelete<IssmDouble>(y_zerolevelset); 4658 xDelete<IssmDouble>(xp); 4659 xDelete<IssmDouble>(yp); 4660 xDelete<IssmDouble>(x); 4661 xDelete<IssmDouble>(y); 4662 xDelete<IssmDouble>(z); 4663 delete vx_zerolevelset; 4664 delete vy_zerolevelset; 4418 4665 } 4419 4666 /*}}}*/ … … 4480 4727 xDelete<IssmDouble>(deviatorictensorerror); 4481 4728 xDelete<IssmDouble>(thicknesserror); 4482 4483 4729 } 4484 4730 /*}}}*/ … … 4512 4758 /*Create initial mesh (coarse mesh) in neopz data structure*/ 4513 4759 /*Just CPU #0 should keep AMR object*/ 4514 this->SetRefPatterns(); 4760 /*Initialize refinement pattern*/ 4761 this->SetRefPatterns(); 4515 4762 if(my_rank==0){ 4516 bool ismisomip = false; 4517 if(ismisomip){//itapopo 4518 TPZFileStream fstr; 4519 std::stringstream ss; 4520 int frictionlaw; 4521 4522 this->parameters->FindParam(&frictionlaw,FrictionLawEnum); 4523 4524 ss << levelmax; 4525 if(frictionlaw==1){ 4526 ss << "_viscous/amr.txt"; 4527 }else if(frictionlaw==7){ 4528 ss << "_tsai/amr.txt"; 4529 }else{ 4530 _error_("friction law not supported here."); 4531 } 4532 4533 std::string AMRfile = "/home/santos/Misomip2/L" + ss.str(); 4534 fstr.OpenRead(AMRfile.c_str()); 4535 4536 TPZSaveable *sv = TPZSaveable::Restore(fstr,0); 4537 this->amr = dynamic_cast<AdaptiveMeshRefinement*>(sv); 4538 } 4539 else{ 4540 this->amr = new AdaptiveMeshRefinement(); 4541 //this->amr->SetLevelMax(levelmax); //Set max level of refinement 4542 //this->amr->SetRegions(regionlevel1,regionlevelmax); 4543 this->amr->CreateInitialMesh(numberofvertices,numberofelements,elementswidth,x,y,elements); 4544 } 4763 this->amr = new AdaptiveMeshRefinement(); 4764 this->amr->CreateInitialMesh(numberofvertices,numberofelements,elementswidth,x,y,elements); 4545 4765 this->amr->SetLevelMax(levelmax); //Set max level of refinement 4546 4766 this->amr->SetRegions(regionlevel1,regionlevelmax); … … 4559 4779 gRefDBase.InitializeUniformRefPattern(ETriangle); 4560 4780 4561 4781 /*Insert specifics patterns to ISSM core*/ 4562 4782 std::string filepath = REFPATTERNDIR; 4563 4783 std::string filename1 = filepath + "/2D_Triang_Rib_3.rpt"; 4564 4784 std::string filename2 = filepath + "/2D_Triang_Rib_4.rpt"; 4565 std::string filename3 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_4.rpt"; 4566 std::string filename4 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_4_permuted.rpt"; 4567 std::string filename5 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_5.rpt"; 4568 std::string filename6 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_5_permuted.rpt"; 4569 std::string filename7 = filepath + "/2D_Triang_Rib_5.rpt"; 4785 std::string filename3 = filepath + "/2D_Triang_Rib_5.rpt"; 4786 std::string filename4 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_4.rpt"; 4787 std::string filename5 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_4_permuted.rpt"; 4788 std::string filename6 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_5.rpt"; 4789 std::string filename7 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_3_5_permuted.rpt"; 4790 std::string filename8 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_4_5.rpt"; 4791 std::string filename9 = filepath + "/2D_Triang_Rib_OnlyTriang_Side_4_5_permuted.rpt"; 4570 4792 4571 4793 TPZAutoPointer<TPZRefPattern> refpat1 = new TPZRefPattern(filename1); … … 4576 4798 TPZAutoPointer<TPZRefPattern> refpat6 = new TPZRefPattern(filename6); 4577 4799 TPZAutoPointer<TPZRefPattern> refpat7 = new TPZRefPattern(filename7); 4800 TPZAutoPointer<TPZRefPattern> refpat8 = new TPZRefPattern(filename8); 4801 TPZAutoPointer<TPZRefPattern> refpat9 = new TPZRefPattern(filename9); 4578 4802 4579 4803 if(!gRefDBase.FindRefPattern(refpat1)) gRefDBase.InsertRefPattern(refpat1); … … 4584 4808 if(!gRefDBase.FindRefPattern(refpat6)) gRefDBase.InsertRefPattern(refpat6); 4585 4809 if(!gRefDBase.FindRefPattern(refpat7)) gRefDBase.InsertRefPattern(refpat7); 4810 if(!gRefDBase.FindRefPattern(refpat8)) gRefDBase.InsertRefPattern(refpat8); 4811 if(!gRefDBase.FindRefPattern(refpat9)) gRefDBase.InsertRefPattern(refpat9); 4586 4812 } 4587 4813 /*}}}*/ -
issm/trunk-jpl/src/c/classes/FemModel.h
r21908 r21919 177 177 void ThicknessZZErrorEstimator(IssmDouble** pelementerror); 178 178 void MeanGroundedIceLevelSet(IssmDouble** pmasklevelset); 179 void GetElementCenterCoordinates(IssmDouble** pxc,IssmDouble** pyc); 179 180 #endif 180 181 … … 182 183 void ReMeshBamg(int* pnewnumberofvertices,int* pnewnumberofelements,IssmDouble** pnewx,IssmDouble** pnewy,IssmDouble** pnewz,int** pnewelementslist); 183 184 void InitializeAdaptiveRefinementBamg(void); 185 void GethmaxVerticesFromZeroLevelSetDistance(IssmDouble* hmaxvertices,int levelset_type); 186 void GethmaxVerticesFromEstimators(IssmDouble* hmaxvertices,int errorestimator_type); 187 void GetVerticeDistanceToZeroLevelSet(IssmDouble** pverticedistance,int leveset_type); 184 188 #endif 185 189 186 190 #if defined(_HAVE_NEOPZ_) && !defined(_HAVE_ADOLC_) 187 /*Adaptive mesh refinement methods*/188 191 void ReMeshNeopz(int* pnewnumberofvertices,int* pnewnumberofelements,IssmDouble** pnewx,IssmDouble** pnewy,IssmDouble** pnewz,int** pnewelementslist); 189 192 void InitializeAdaptiveRefinementNeopz(void); -
issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp
r21830 r21919 114 114 parameters->AddObject(iomodel->CopyConstantObject("md.amr.hmax",AmrHmaxEnum)); 115 115 parameters->AddObject(iomodel->CopyConstantObject("md.amr.err",AmrErrEnum)); 116 parameters->AddObject(iomodel->CopyConstantObject("md.amr.keepmetric",AmrKeepMetricEnum)); 117 parameters->AddObject(iomodel->CopyConstantObject("md.amr.groundingline_resolution",AmrGroundingLineResolutionEnum)); 118 parameters->AddObject(iomodel->CopyConstantObject("md.amr.groundingline_distance",AmrGroundingLineDistanceEnum)); 119 parameters->AddObject(iomodel->CopyConstantObject("md.amr.icefront_resolution",AmrIceFrontResolutionEnum)); 120 parameters->AddObject(iomodel->CopyConstantObject("md.amr.icefront_distance",AmrIceFrontDistanceEnum)); 121 parameters->AddObject(iomodel->CopyConstantObject("md.amr.thicknesserror_resolution",AmrThicknessErrorResolutionEnum)); 122 parameters->AddObject(iomodel->CopyConstantObject("md.amr.thicknesserror_threshold",AmrThicknessErrorThresholdEnum)); 123 parameters->AddObject(iomodel->CopyConstantObject("md.amr.deviatoricerror_resolution",AmrDeviatoricErrorResolutionEnum)); 124 parameters->AddObject(iomodel->CopyConstantObject("md.amr.deviatoricerror_threshold",AmrDeviatoricErrorThresholdEnum)); 116 125 /*Convert fieldname to enum and put it in params*/ 117 126 iomodel->FindConstant(&fieldname,"md.amr.fieldname"); -
issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
r21908 r21919 849 849 AmrFieldEnum, 850 850 AmrErrEnum, 851 AmrKeepMetricEnum, 852 AmrGroundingLineResolutionEnum, 853 AmrGroundingLineDistanceEnum, 854 AmrIceFrontResolutionEnum, 855 AmrIceFrontDistanceEnum, 856 AmrThicknessErrorResolutionEnum, 857 AmrThicknessErrorThresholdEnum, 858 AmrDeviatoricErrorResolutionEnum, 859 AmrDeviatoricErrorThresholdEnum, 851 860 DeviatoricStressErrorEstimatorEnum, 852 861 ThicknessErrorEstimatorEnum, -
issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
r21908 r21919 825 825 case AmrFieldEnum : return "AmrField"; 826 826 case AmrErrEnum : return "AmrErr"; 827 case AmrKeepMetricEnum : return "AmrKeepMetric"; 828 case AmrGroundingLineResolutionEnum : return "AmrGroundingLineResolution"; 829 case AmrGroundingLineDistanceEnum : return "AmrGroundingLineDistance"; 830 case AmrIceFrontResolutionEnum : return "AmrIceFrontResolution"; 831 case AmrIceFrontDistanceEnum : return "AmrIceFrontDistance"; 832 case AmrThicknessErrorResolutionEnum : return "AmrThicknessErrorResolution"; 833 case AmrThicknessErrorThresholdEnum : return "AmrThicknessErrorThreshold"; 834 case AmrDeviatoricErrorResolutionEnum : return "AmrDeviatoricErrorResolution"; 835 case AmrDeviatoricErrorThresholdEnum : return "AmrDeviatoricErrorThreshold"; 827 836 case DeviatoricStressErrorEstimatorEnum : return "DeviatoricStressErrorEstimator"; 828 837 case ThicknessErrorEstimatorEnum : return "ThicknessErrorEstimator"; -
issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
r21908 r21919 843 843 else if (strcmp(name,"AmrField")==0) return AmrFieldEnum; 844 844 else if (strcmp(name,"AmrErr")==0) return AmrErrEnum; 845 else if (strcmp(name,"AmrKeepMetric")==0) return AmrKeepMetricEnum; 846 else if (strcmp(name,"AmrGroundingLineResolution")==0) return AmrGroundingLineResolutionEnum; 847 else if (strcmp(name,"AmrGroundingLineDistance")==0) return AmrGroundingLineDistanceEnum; 848 else if (strcmp(name,"AmrIceFrontResolution")==0) return AmrIceFrontResolutionEnum; 849 else if (strcmp(name,"AmrIceFrontDistance")==0) return AmrIceFrontDistanceEnum; 850 else if (strcmp(name,"AmrThicknessErrorResolution")==0) return AmrThicknessErrorResolutionEnum; 851 else if (strcmp(name,"AmrThicknessErrorThreshold")==0) return AmrThicknessErrorThresholdEnum; 852 else if (strcmp(name,"AmrDeviatoricErrorResolution")==0) return AmrDeviatoricErrorResolutionEnum; 853 else if (strcmp(name,"AmrDeviatoricErrorThreshold")==0) return AmrDeviatoricErrorThresholdEnum; 845 854 else if (strcmp(name,"DeviatoricStressErrorEstimator")==0) return DeviatoricStressErrorEstimatorEnum; 846 855 else if (strcmp(name,"ThicknessErrorEstimator")==0) return ThicknessErrorEstimatorEnum; … … 866 875 else if (strcmp(name,"ExternalResult")==0) return ExternalResultEnum; 867 876 else if (strcmp(name,"FileParam")==0) return FileParamEnum; 868 else if (strcmp(name,"Input")==0) return InputEnum; 877 else stage=8; 878 } 879 if(stage==8){ 880 if (strcmp(name,"Input")==0) return InputEnum; 869 881 else if (strcmp(name,"IntInput")==0) return IntInputEnum; 870 882 else if (strcmp(name,"IntParam")==0) return IntParamEnum; … … 875 887 else if (strcmp(name,"Matenhancedice")==0) return MatenhancediceEnum; 876 888 else if (strcmp(name,"Matestar")==0) return MatestarEnum; 877 else stage=8; 878 } 879 if(stage==8){ 880 if (strcmp(name,"Matpar")==0) return MatparEnum; 889 else if (strcmp(name,"Matpar")==0) return MatparEnum; 881 890 else if (strcmp(name,"Node")==0) return NodeEnum; 882 891 else if (strcmp(name,"Numericalflux")==0) return NumericalfluxEnum; … … 989 998 else if (strcmp(name,"MaxVel")==0) return MaxVelEnum; 990 999 else if (strcmp(name,"MinVx")==0) return MinVxEnum; 991 else if (strcmp(name,"MaxVx")==0) return MaxVxEnum; 1000 else stage=9; 1001 } 1002 if(stage==9){ 1003 if (strcmp(name,"MaxVx")==0) return MaxVxEnum; 992 1004 else if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum; 993 1005 else if (strcmp(name,"MinVy")==0) return MinVyEnum; … … 998 1010 else if (strcmp(name,"MaxAbsVz")==0) return MaxAbsVzEnum; 999 1011 else if (strcmp(name,"FloatingArea")==0) return FloatingAreaEnum; 1000 else stage=9; 1001 } 1002 if(stage==9){ 1003 if (strcmp(name,"GroundedArea")==0) return GroundedAreaEnum; 1012 else if (strcmp(name,"GroundedArea")==0) return GroundedAreaEnum; 1004 1013 else if (strcmp(name,"IceMass")==0) return IceMassEnum; 1005 1014 else if (strcmp(name,"IceVolume")==0) return IceVolumeEnum; -
issm/trunk-jpl/src/m/classes/amrbamg.m
r21802 r21919 10 10 fieldname = ''; 11 11 err = 0.; 12 keepmetric = 0; 13 groundingline_resolution = 0.; 14 groundingline_distance = 0.; 15 icefront_resolution = 0.; 16 icefront_distance = 0.; 17 thicknesserror_resolution = 0.; 18 thicknesserror_threshold = 0.; 19 deviatoricerror_resolution = 0.; 20 deviatoricerror_threshold = 0.; 12 21 end 13 22 methods … … 30 39 self.err=3.; 31 40 41 %keep metric? 42 self.keepmetric=1; 43 44 %other criterias 45 self.groundingline_resolution=500.; 46 self.groundingline_distance=0.; 47 self.icefront_resolution=500.; 48 self.icefront_distance=0.; 49 self.thicknesserror_resolution=500.; 50 self.thicknesserror_threshold=0.; 51 self.deviatoricerror_resolution=500.; 52 self.deviatoricerror_threshold=0.; 53 32 54 end % }}} 33 55 function md = checkconsistency(self,md,solution,analyses) % {{{ … … 36 58 md = checkfield(md,'fieldname','amr.hmin','numel',[1],'>',0,'<',self.hmax,'NaN',1); 37 59 %md = checkfield(md,'fieldname','amr.fieldname','string',[1]); 60 md = checkfield(md,'fieldname','amr.keepmetric','numel',[1],'>=',0,'<=',1,'NaN',1); 61 md = checkfield(md,'fieldname','amr.groundingline_resolution','numel',[1],'>',0,'<',self.hmax,'NaN',1); 62 md = checkfield(md,'fieldname','amr.groundingline_distance','numel',[1],'>=',0,'NaN',1,'Inf',1); 63 md = checkfield(md,'fieldname','amr.icefront_resolution','numel',[1],'>',0,'<',self.hmax,'NaN',1); 64 md = checkfield(md,'fieldname','amr.icefront_distance','numel',[1],'>=',0,'NaN',1,'Inf',1); 65 md = checkfield(md,'fieldname','amr.thicknesserror_resolution','numel',[1],'>',0,'<',self.hmax,'NaN',1); 66 md = checkfield(md,'fieldname','amr.thicknesserror_threshold','numel',[1],'>=',0,'<=',1,'NaN',1); 67 md = checkfield(md,'fieldname','amr.deviatoricerror_resolution','numel',[1],'>',0,'<',self.hmax,'NaN',1); 68 md = checkfield(md,'fieldname','amr.deviatoricerror_threshold','numel',[1],'>=',0,'<=',1,'NaN',1); 38 69 end % }}} 39 70 function disp(self) % {{{ … … 43 74 fielddisplay(self,'hmax',['maximum element length']); 44 75 fielddisplay(self,'fieldname',['name of input that will be used to compute the metric (should be an input of FemModel)']); 76 fielddisplay(self,'keepmetric',['indicates whether the metric should be kept every remeshing time']); 77 fielddisplay(self,'groundingline_resolution',['element length near the grounding line']); 78 fielddisplay(self,'groundingline_distance',['distance around the grounding line which elements will be refined']); 79 fielddisplay(self,'icefront_resolution',['element length near the ice front']); 80 fielddisplay(self,'icefront_distance',['distance around the ice front which elements will be refined']); 81 fielddisplay(self,'thicknesserror_resolution',['element length when thickness error estimator is used']); 82 fielddisplay(self,'thicknesserror_threshold',['maximum threshold thickness error permitted']); 83 fielddisplay(self,'deviatoricerror_resolution',['element length when deviatoric stress error estimator is used']); 84 fielddisplay(self,'deviatoricerror_threshold',['maximum threshold deviatoricstress error permitted']); 45 85 46 86 end % }}} … … 51 91 WriteData(fid,prefix,'object',self,'class','amr','fieldname','hmax','format','Double'); 52 92 WriteData(fid,prefix,'object',self,'class','amr','fieldname','fieldname','format','String'); 53 WriteData(fid,prefix,'object',self,'class','err','fieldname','err','format','Double'); 93 WriteData(fid,prefix,'object',self,'class','amr','fieldname','err','format','Double'); 94 WriteData(fid,prefix,'object',self,'class','amr','fieldname','keepmetric','format','Integer'); 95 WriteData(fid,prefix,'object',self,'class','amr','fieldname','groundingline_resolution','format','Double'); 96 WriteData(fid,prefix,'object',self,'class','amr','fieldname','groundingline_distance','format','Double'); 97 WriteData(fid,prefix,'object',self,'class','amr','fieldname','icefront_resolution','format','Double'); 98 WriteData(fid,prefix,'object',self,'class','amr','fieldname','icefront_distance','format','Double'); 99 WriteData(fid,prefix,'object',self,'class','amr','fieldname','thicknesserror_resolution','format','Double'); 100 WriteData(fid,prefix,'object',self,'class','amr','fieldname','thicknesserror_threshold','format','Double'); 101 WriteData(fid,prefix,'object',self,'class','amr','fieldname','deviatoricerror_resolution','format','Double'); 102 WriteData(fid,prefix,'object',self,'class','amr','fieldname','deviatoricerror_threshold','format','Double'); 54 103 55 104 end % }}} -
issm/trunk-jpl/src/m/mesh/bamg.m
r21864 r21919 25 25 % - maxsubdiv : maximum subdivision of exisiting elements (default is 10) 26 26 % - metric : matrix (numberofnodes x 3) used as a metric 27 % - Metrictype : 1-> absolute error c/(err coeff^2) * Abs(H) (default)28 % 2-> relative error c/(err coeff^2) * Abs(H)/max(s,cutoff*max(s))29 % 3-> rescaled absolute error c/(err coeff^2) * Abs(H)/(smax-smin)27 % - Metrictype : 0 -> absolute error c/(err coeff^2) * Abs(H) (default) 28 % 1 -> relative error c/(err coeff^2) * Abs(H)/max(s,cutoff*max(s)) 29 % 2 -> rescaled absolute error c/(err coeff^2) * Abs(H)/(smax-smin) 30 30 % - nbjacoby : correction used by Hessiantype=1 (default is 1) 31 31 % - nbsmooth : number of metric smoothing procedure (default is 3)
Note:
See TracChangeset
for help on using the changeset viewer.
