Changeset 21806


Ignore:
Timestamp:
07/18/17 08:05:13 (8 years ago)
Author:
tsantos
Message:

CHG: added error estimator for thickness, adjuted AMR with NeoPZ. Working with unrefinement and delete element process.

Location:
issm/trunk-jpl/src
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp

    r21762 r21806  
    1010
    1111#include "./AdaptiveMeshRefinement.h"
     12#include "TPZVTKGeoMesh.h"
     13#include "../shared/shared.h"
     14#include "pzgeotriangle.h"
     15#include "pzreftriangle.h"
     16using namespace pzgeom;
    1217
    1318/*Constructor, copy, clean up and destructor*/
    … …  
    1722/*}}}*/
    1823AdaptiveMeshRefinement::AdaptiveMeshRefinement(const AdaptiveMeshRefinement &cp){/*{{{*/
    19    this->Initialize();
     24        this->Initialize();
    2025        this->operator =(cp);
    2126}
    … …  
    2530        /*Clean all attributes*/
    2631        this->CleanUp();
    27 
    2832        /*Copy all data*/
    29         this->fathermesh     = new TPZGeoMesh(*cp.fathermesh);
    30         this->previousmesh   = new TPZGeoMesh(*cp.previousmesh);
    31         this->levelmax       = cp.levelmax;
    32         this->elementswidth  = cp.elementswidth;
    33         this->regionlevel1   = cp.regionlevel1;
    34         this->regionlevelmax = cp.regionlevelmax;
     33        this->fathermesh                        = new TPZGeoMesh(*cp.fathermesh);
     34        this->currentmesh                       = new TPZGeoMesh(*cp.currentmesh);
     35        this->levelmax                          = cp.levelmax;
     36        this->elementswidth             = cp.elementswidth;
     37        this->regionlevel1              = cp.regionlevel1;
     38        this->regionlevelmax            = cp.regionlevelmax;
     39        this->sid2index.clear();
     40        this->sid2index.resize(cp.sid2index.size());
     41        for(int i=0;i<cp.sid2index.size();i++) this->sid2index[i]=cp.sid2index[i];
     42       
    3543        return *this;
    3644
    … …  
    3947AdaptiveMeshRefinement::~AdaptiveMeshRefinement(){/*{{{*/
    4048       
    41         //bool ismismip = false;
    42         //if(ismismip){//itapopo
    43         //      TPZFileStream fstr;
    44         //      std::stringstream ss;
     49        bool ismismip = false;
     50        if(ismismip){//itapopo
     51                TPZFileStream fstr;
     52                std::stringstream ss;
    4553           
    46         //      ss << this->levelmax;
    47         //      std::string AMRfile     = "/home/santos/Misomip2/L" + ss.str() + "_tsai/amr.txt";
    48        
    49         //      fstr.OpenWrite(AMRfile.c_str());
    50         //      int withclassid = 1;
    51         //      this->Write(fstr,withclassid);
    52         //}
     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        }
    5361        this->CleanUp();
    5462        gRefDBase.clear();
    … …  
    5765void AdaptiveMeshRefinement::CleanUp(){/*{{{*/
    5866
    59     /*Verify and delete all data*/
     67        /*Verify and delete all data*/
    6068        if(this->fathermesh)    delete this->fathermesh;
    61    if(this->previousmesh)  delete this->previousmesh;
    62         this->levelmax                  = -1;
    63         this->elementswidth  = -1;
    64         this->regionlevel1      = -1;
    65         this->regionlevelmax = -1;
     69        if(this->currentmesh)   delete this->currentmesh;
     70        this->levelmax                          = -1;
     71        this->elementswidth             = -1;
     72        this->regionlevel1              = -1;
     73        this->regionlevelmax            = -1;
     74        this->sid2index.clear();
    6675}
    6776/*}}}*/
    … …  
    6978
    7079        /*Set pointers to NULL*/
    71         this->fathermesh                = NULL;
    72         this->previousmesh      = NULL;
    73         this->levelmax                  = -1;
    74         this->elementswidth     = -1;
    75         this->regionlevel1      = -1;
    76         this->regionlevelmax = -1;
     80        this->fathermesh                        = NULL;
     81        this->currentmesh                       = NULL;
     82        this->levelmax                          = -1;
     83        this->elementswidth             = -1;
     84        this->regionlevel1              = -1;
     85        this->regionlevelmax            = -1;
     86        this->sid2index.clear();
    7787}
    7888/*}}}*/
    … …  
    8191}
    8292/*}}}*/
    83 void AdaptiveMeshRefinement::Read(TPZStream &buf, void *context){/*{{{*/
    84 
    85     try
    86     {
    87         /* Read the id context*/
    88         TPZSaveable::Read(buf,context);
    89 
    90         /* Read class id*/
    91         int classid;
    92         buf.Read(&classid,1);
    93        
    94         /* Verify the class id*/
    95         if (classid != this->ClassId() )
    96         {
    97             std::cout << "Error in restoring AdaptiveMeshRefinement!\n";
    98             std::cout.flush();
    99             DebugStop();
    100         }
    101        
    102         /* Read simple attributes */
    103         buf.Read(&this->levelmax,1);
    104         buf.Read(&this->elementswidth,1);
    105         buf.Read(&this->regionlevel1,1);
    106         buf.Read(&this->regionlevelmax,1);
    107        
    108                 /* Read geometric mesh*/
    109         TPZSaveable *sv1 = TPZSaveable::Restore(buf,0);
    110         this->fathermesh = dynamic_cast<TPZGeoMesh*>(sv1);
    111        
    112         TPZSaveable *sv2 = TPZSaveable::Restore(buf,0);
    113         this->previousmesh = dynamic_cast<TPZGeoMesh*>(sv2);
     93void 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/*}}}*/
     132template class TPZRestoreClass<AdaptiveMeshRefinement,13829430>;/*{{{*/
     133/*}}}*/
     134void 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());
    114156    }
    115157    catch(const std::exception& e)
    116158    {
    117         std::cout << "Exception catched! " << e.what() << std::endl;
    118         std::cout.flush();
    119         DebugStop();
     159                _error_("AdaptiveMeshRefinement::Write: Exception catched!\n");
    120160    }
    121161}
    122162/*}}}*/
    123 template class TPZRestoreClass<AdaptiveMeshRefinement,13829430>;/*{{{*/
    124 /*}}}*/
    125 void AdaptiveMeshRefinement::Write(TPZStream &buf, int withclassid){/*{{{*/
    126    
    127     try
    128     {
    129         /* Write context (this class) class ID*/
    130         TPZSaveable::Write(buf,withclassid);
    131 
    132         /* Write this class id*/
    133         int classid = ClassId();
    134         buf.Write(&classid,1);
    135 
    136         /* Write simple attributes */
    137         buf.Write(&this->levelmax,1);
    138         buf.Write(&this->elementswidth,1);
    139         buf.Write(&this->regionlevel1,1);
    140         buf.Write(&this->regionlevelmax,1);
    141                        
    142         /* Write the geometric mesh*/
    143         this->fathermesh->Write(buf, this->ClassId());
    144         this->previousmesh->Write(buf, this->ClassId());
    145     }
    146     catch(const std::exception& e)
    147     {
    148         std::cout << "Exception catched! " << e.what() << std::endl;
    149         std::cout.flush();
    150         DebugStop();
    151     }
    152 }
    153 /*}}}*/
    154163
    155164/*Mesh refinement methods*/
    156 #include "TPZVTKGeoMesh.h" //itapopo
    157 #include "../shared/shared.h" //itapopo
    158 void AdaptiveMeshRefinement::ExecuteRefinement(int &type_process,double *vx, double *vy, double *masklevelset, int &nvertices, int &nelements, int &nsegments, double** px, double** py, double** pz, int** pelements, int** psegments){/*{{{*/
    159 
    160         /*IMPORTANT! pelements (and psegments) are in Matlab indexing*/
     165void AdaptiveMeshRefinement::Execute(bool &amr_verbose,
     166                                                                                                int &numberofelements,
     167                                                                                                double* partiallyfloatedelements,
     168                                                                                                double *masklevelset,
     169                                                                                                double* deviatorictensorerror,
     170                                                                                                double* thicknesserror,
     171                                                                                                int &newnumberofvertices,int &newnumberofelements,double** x,double** y,int** elementslist){/*{{{*/
     172
     173        /*IMPORTANT! newelements are in Matlab indexing*/
    161174        /*NEOPZ works only in C indexing*/
    162 
    163     _assert_(this->fathermesh);
    164     _assert_(this->previousmesh);
    165    
    166     /*Calculate the position of the grounding line using previous mesh*/
    167     std::vector<TPZVec<REAL> > GLvec;
    168     this->CalcGroundingLinePosition(masklevelset, GLvec);
    169    
    170    // std::ofstream file1("/home/santos/mesh0.vtk");
    171    // TPZVTKGeoMesh::PrintGMeshVTK(this->fathermesh,file1 );
    172    
    173     /*run refinement or unrefinement process*/
    174     TPZGeoMesh *newmesh;
    175     switch (type_process) {
    176         case 0: newmesh = this->previousmesh; break;                    // refine previous mesh
    177         case 1: newmesh = new TPZGeoMesh(*this->fathermesh); break;     // refine mesh 0 (unrefine process)
    178         default: DebugStop(); break;//itapopo verificar se irá usar _assert_
    179     }
    180    
    181     this->RefinementProcess(newmesh,GLvec);
    182        
    183     //std::ofstream file2("/home/santos/mesh1.vtk");
    184     //TPZVTKGeoMesh::PrintGMeshVTK(this->previousmesh,file2 );
    185    
    186     /*Set new mesh pointer. Previous mesh just have uniform elements*/
    187     if(type_process==1){
    188         if(this->previousmesh) delete this->previousmesh;
    189         this->previousmesh = newmesh;
    190     }
    191    
    192     /*Refine elements to avoid hanging nodes*/
    193         //TPZGeoMesh *nohangingnodesmesh = new TPZGeoMesh(*newmesh);//itapopo testando, este era o original
    194    TPZGeoMesh *nohangingnodesmesh = this->CreateRefPatternMesh(newmesh);//itapopo testando, este eh novo metodo
    195    
    196     //std::ofstream file3("/home/santos/mesh2.vtk");
    197     //TPZVTKGeoMesh::PrintGMeshVTK(this->previousmesh,file3);
    198    
    199     this->RefineMeshToAvoidHangingNodes(nohangingnodesmesh);
    200    
    201          //std::ofstream file4("/home/santos/mesh3.vtk");
    202     //TPZVTKGeoMesh::PrintGMeshVTK(nohangingnodesmesh,file4);
    203    
    204     /*Get new geometric mesh in ISSM data structure*/
    205     this->GetMesh(nohangingnodesmesh,nvertices,nelements,nsegments,px,py,pz,pelements,psegments);
    206          
    207     /*Verify the new geometry*/
    208     this->CheckMesh(nvertices,nelements,nsegments,this->elementswidth,px,py,pz,pelements,psegments);
    209 
    210          _printf_("\trefinement process done!\n\n");
    211 
    212     delete nohangingnodesmesh;
    213 }
    214 /*}}}*/
    215 void AdaptiveMeshRefinement::RefinementProcess(TPZGeoMesh *gmesh,std::vector<TPZVec<REAL> > &GLvec){/*{{{*/
    216    
    217     /*Refine mesh levelmax times*/
    218    _printf_("\n\trefinement process (level max = " << this->levelmax << ")\n");
    219         _printf_("\tprogress:  ");
    220         for(int hlevel=1;hlevel<=this->levelmax;hlevel++){
    221        
    222         /*Set elements to be refined using some criteria*/
    223         std::vector<int> ElemVec; //elements without children
    224         this->SetElementsToRefine(gmesh,GLvec,hlevel,ElemVec);
    225        
    226         /*Refine the mesh*/
    227         this->RefineMesh(gmesh, ElemVec);
    228                  
    229                   _printf_("*  ");
    230     }
    231     _printf_("\n");
    232 }
    233 /*}}}*/
    234 void AdaptiveMeshRefinement::RefineMesh(TPZGeoMesh *gmesh, std::vector<int> &ElemVec){/*{{{*/
    235 
    236         /*Refine elements in ElemVec: uniform pattern refinement*/
    237         for(int i = 0; i < ElemVec.size(); i++){
    238                
    239         /*Get geometric element and verify if it has already been refined*/
    240         int index = ElemVec[i];
    241         TPZGeoEl * geoel = gmesh->Element(index);
    242         if(geoel->HasSubElement()) DebugStop();                              //itapopo _assert_(!geoel->HasSubElement());
    243         if(geoel->MaterialId() != this->GetElemMaterialID()) DebugStop();   //itapopo verificar se usará _assert_
    244        
    245         /*Divide geoel*/
    246         TPZVec<TPZGeoEl *> Sons;
    247                   geoel->Divide(Sons);
    248        
    249         /*If a 1D segment is neighbor, it must be divided too*/
    250         if(this->elementswidth != 3) DebugStop(); //itapopo verificar o segment para malha 3D
    251        
    252         std::vector<int> sides(3);
    253         sides[0] = 3; sides[1] = 4; sides[2] = 5;
    254         for(int j = 0; j < sides.size(); j++ ){
    255            
    256             TPZGeoElSide Neighbour = geoel->Neighbour(sides[j]);
    257            
    258             if( Neighbour.Element()->MaterialId() == this->GetBoundaryMaterialID() && !Neighbour.Element()->HasSubElement() ){
    259                 TPZVec<TPZGeoEl *> pv2;
    260                 Neighbour.Element()->Divide(pv2);
    261             }
    262         }
    263         }
    264    
    265     gmesh->BuildConnectivity();
    266 
     175        if(!this->fathermesh || !this->currentmesh) _error_("Impossible to execute refinement: fathermesh or currentmesh is NULL!\n");
     176        if(numberofelements!=this->sid2index.size()) _error_("Impossible to execute refinement: sid2index.size is not equal to numberofelements!\n");
     177
     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/*}}}*/
     189void 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;
     196        double mean_mask                = 0;
     197        double mean_tauerror = 0;
     198        double mean_Herror      = 0;
     199        double group_error      = 0;
     200   
     201        /*Calculate mean values*/
     202        for(int i=0;i<this->sid2index.size();i++){
     203                mean_mask               += masklevelset[i];
     204                mean_tauerror   += deviatorictensorerror[i];
     205                mean_Herror             += thicknesserror[i];
     206        }
     207        mean_mask               /= this->sid2index.size();
     208        mean_tauerror   /= this->sid2index.size();
     209        mean_Herror             /= this->sid2index.size();
     210
     211        if(amr_verbose) _printf_("\t\tuniform refinement...\n");
     212        for(int i=0;i<this->sid2index.size();i++){
     213                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){
     218                        TPZVec<TPZGeoEl *> sons;
     219                        if(geoel->Level()<this->levelmax) geoel->Divide(sons);
     220                }
     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();
     227                        }
     228                }
     229        }
     230        this->currentmesh->BuildConnectivity();
     231       
     232        if(amr_verbose) _printf_("\t\trefine to avoid hanging nodes...\n");
     233        this->RefineMeshToAvoidHangingNodes(this->currentmesh);
     234        this->currentmesh->BuildConnectivity();
     235       
     236                //nohangingnodesmesh = this->CreateRefPatternMesh(newmesh); itapopo tentar otimizar
     237       
     238        if(amr_verbose) _printf_("\trefinement process done!\n");
     239}
     240/*}}}*/
     241void AdaptiveMeshRefinement::RefineMesh(TPZGeoMesh *gmesh,std::vector<int> &elements){/*{{{*/
     242
     243        /*Refine elements: uniform pattern refinement*/
     244        for(int i=0;i<elements.size();i++){
     245                /*Get geometric element and verify if it has already been refined*/
     246                int index = elements[i];
     247                TPZGeoEl * geoel = gmesh->Element(index);
     248                if(geoel->HasSubElement()) _error_("Impossible to refine: geoel (index) " << index << " has subelements!\n");
     249                if(geoel->MaterialId()!=this->GetElemMaterialID()) _error_("Impossible to refine: geoel->MaterialId is not GetElemMaterialID!\n");
     250                /*Divide geoel*/
     251                TPZVec<TPZGeoEl *> Sons;
     252                geoel->Divide(Sons);
     253        }
     254        gmesh->BuildConnectivity();
    267255}
    268256/*}}}*/
    269257void AdaptiveMeshRefinement::RefineMeshToAvoidHangingNodes(TPZGeoMesh *gmesh){/*{{{*/
    270258   
    271          _printf_("\trefine to avoid hanging nodes...\n");
    272     /*Refine elements to avoid hanging nodes: non-uniform refinement*/
    273          const int NElem = gmesh->NElements();
    274     for(int i = 0; i < NElem; i++){
    275        
    276         /*Get geometric element and verify if it has already been refined. Geoel may not have been previously refined*/
    277         TPZGeoEl * geoel = gmesh->Element(i);
    278         if(!geoel) continue;
    279         if(geoel->HasSubElement()) continue;
    280         if(geoel->MaterialId() != this->GetElemMaterialID()) continue;
    281        
    282         /*Get the refinement pattern for this element and refine it*/
    283         TPZAutoPointer<TPZRefPattern> refp = TPZRefPatternTools::PerfectMatchRefPattern(geoel);
    284         if(refp){
    285             TPZVec<TPZGeoEl *> Sons;
    286             geoel->SetRefPattern(refp);
    287             geoel->Divide(Sons);
    288         }
    289        
    290     }
    291    
    292     gmesh->BuildConnectivity();
    293    
    294 }
    295 /*}}}*/
    296 void AdaptiveMeshRefinement::GetMesh(TPZGeoMesh *gmesh,int &nvertices,int &nelements,int &nsegments,double** px,double** py,double** pz, int** pelements, int** psegments){/*{{{*/
    297 
    298         /*IMPORTANT! pelements (and psegments) are in Matlab indexing*/
    299         /*NEOPZ works only in C indexing*/
    300 
    301         /* vertices */
    302     int ntotalvertices = gmesh->NNodes();//total
    303    
    304     /* mesh coords */
    305          double* newmeshX = xNew<IssmDouble>(ntotalvertices);
    306     double* newmeshY = xNew<IssmDouble>(ntotalvertices);
    307     double* newmeshZ = xNew<IssmDouble>(ntotalvertices);
    308        
    309    /* getting mesh coords */
    310     for(int i = 0; i < ntotalvertices; i++ ){
    311         TPZVec<REAL> coords(3,0.);
    312         gmesh->NodeVec()[i].GetCoordinates(coords);
    313         newmeshX[i] = coords[0];
    314         newmeshY[i] = coords[1];
    315         newmeshZ[i] = coords[2];
    316     }
    317    
    318         /* elements */
    319     std::vector<TPZGeoEl*> GeoVec; GeoVec.clear();
    320     for(int i = 0; i < gmesh->NElements(); i++){
    321         if( gmesh->ElementVec()[i]->HasSubElement() ) continue;
    322         if( gmesh->ElementVec()[i]->MaterialId() != this->GetElemMaterialID() ) continue;
    323         GeoVec.push_back( gmesh->ElementVec()[i]);
    324     }
    325    
    326     int ntotalelements = (int)GeoVec.size();
    327     int* newelements = xNew<int>(ntotalelements*this->elementswidth);
    328 
    329     if ( !(this->elementswidth == 3) && !(this->elementswidth == 4) && !(this->elementswidth == 6) ) DebugStop();
    330 
    331     for(int i=0;i<GeoVec.size();i++){
    332         for(int j=0;j<this->elementswidth;j++) newelements[i*this->elementswidth+j]=(int)GeoVec[i]->NodeIndex(j)+1;//C to Matlab indexing
    333          }
    334    
    335     /* segments */
    336     std::vector<TPZGeoEl*> SegVec; SegVec.clear();
    337     for(int i = 0; i < gmesh->NElements(); i++){
    338         if( gmesh->ElementVec()[i]->HasSubElement() ) continue;
    339         if( gmesh->ElementVec()[i]->MaterialId() != this->GetBoundaryMaterialID() ) continue;
    340         SegVec.push_back( gmesh->ElementVec()[i]);
    341     }
    342    
    343     int ntotalsegments = (int)SegVec.size();
    344     int *newsegments=NULL;
    345          if(ntotalsegments>0) newsegments=xNew<int>(ntotalsegments*3);
    346    
    347     for(int i=0;i<SegVec.size();i++){
    348        
    349         for(int j=0;j<2;j++) newsegments[i*3+j]=(int)SegVec[i]->NodeIndex(j)+1;//C to Matlab indexing
    350        
    351         int neighborindex = SegVec[i]->Neighbour(2).Element()->Index();
    352         int neighbourid = -1;
    353        
    354         for(int j = 0; j < GeoVec.size(); j++){
    355             if( GeoVec[j]->Index() == neighborindex || GeoVec[j]->FatherIndex() == neighborindex){
    356                 neighbourid = j;
    357                 break;
    358             }
    359         }
    360        
    361         if(neighbourid==-1) DebugStop(); //itapopo talvez passar para _assert_
    362         newsegments[i*3+2] = neighbourid+1;//C to Matlab indexing
    363     }
    364    
    365     //setting outputs
    366     nvertices  = ntotalvertices;
    367     nelements  = ntotalelements;
    368     nsegments  = ntotalsegments;
    369     *px            = newmeshX;
    370     *py            = newmeshY;
    371     *pz            = newmeshZ;
    372     *pelements = newelements;
    373     *psegments = newsegments;
    374    
    375 }
    376 /*}}}*/
    377 void AdaptiveMeshRefinement::CalcGroundingLinePosition(double *masklevelset,std::vector<TPZVec<REAL> > &GLvec){/*{{{*/
    378    
    379     /* Find grounding line using elments center point */
    380     GLvec.clear();
    381     for(int i=0;i<this->previousmesh->NElements();i++){
    382        
    383         if(this->previousmesh->Element(i)->MaterialId()!=this->GetElemMaterialID()) continue;
    384         if(this->previousmesh->Element(i)->HasSubElement()) continue;
    385        
    386         //itapopo apenas malha 2D triangular!
    387         int vertex0 = this->previousmesh->Element(i)->NodeIndex(0);
    388         int vertex1 = this->previousmesh->Element(i)->NodeIndex(1);
    389         int vertex2 = this->previousmesh->Element(i)->NodeIndex(2);
    390        
    391                   //itapopo inserir uma verificação para não acessar fora da memória
    392         double mls0 = masklevelset[vertex0];
    393         double mls1 = masklevelset[vertex1];
    394         double mls2 = masklevelset[vertex2];
    395        
    396         if( mls0*mls1 < 0. || mls1*mls2 < 0. ){
    397             const int side = 6;
    398             TPZVec<double> qsi(2,0.);
    399             TPZVec<double> X(3,0.);
    400             this->previousmesh->Element(i)->CenterPoint(side, qsi);
    401             this->previousmesh->Element(i)->X(qsi, X);
    402             GLvec.push_back(X);
    403         }
    404     }
    405    
    406 //    itapopo apenas para debugar
    407 //    std::ofstream fileGL("/Users/santos/Desktop/gl.nb");
    408 //    fileGL << "ListPlot[{";
    409 //    for(int i = 0; i < GLvec.size(); i++){
    410 //        fileGL << "{" << GLvec[i][0] << "," << GLvec[i][1] << /*"," << 0. << */"}";
    411 //        if(i != GLvec.size()-1) fileGL << ",";
    412 //    }
    413 //    fileGL << "}]";
    414 //    fileGL.flush();
    415 //    fileGL.close();
    416    
    417 }
    418 /*}}}*/
    419 void AdaptiveMeshRefinement::SetElementsToRefine(TPZGeoMesh *gmesh,std::vector<TPZVec<REAL> > &GLvec,int &hlevel,std::vector<int> &ElemVec){/*{{{*/
    420 
    421     if(!gmesh) DebugStop(); //itapopo verificar se usará _assert_
    422 
    423     ElemVec.clear();
    424    
    425     // itapopo inserir modo de encontrar criterio TESTING!!!! Come to false!
    426          if(false) this->TagAllElements(gmesh,ElemVec); //uniform, refine all elements!
    427 
    428     /* Adaptive refinement. This refines some elements following some criteria*/
    429     this->TagElementsNearGroundingLine(gmesh, GLvec, hlevel, ElemVec);
    430 
    431 }
    432 /*}}}*/
    433 void AdaptiveMeshRefinement::TagAllElements(TPZGeoMesh *gmesh,std::vector<int> &ElemVec){/*{{{*/
    434    
     259        /*Refine elements to avoid hanging nodes: non-uniform refinement*/
     260        const int NElem = gmesh->NElements();
     261        for(int i=0;i<NElem;i++){
     262                /*Get geometric element and verify if it has already been refined. Geoel may not have been previously refined*/
     263                TPZGeoEl * geoel=gmesh->Element(i);
     264                if(!geoel) continue;
     265                if(geoel->HasSubElement()) continue;
     266                if(geoel->MaterialId() != this->GetElemMaterialID()) continue;
     267                /*Get the refinement pattern for this element and refine it*/
     268                TPZAutoPointer<TPZRefPattern> refp=TPZRefPatternTools::PerfectMatchRefPattern(geoel);
     269                if(refp){
     270                        TPZVec<TPZGeoEl *> Sons;
     271                        geoel->SetRefPattern(refp);
     272                        geoel->Divide(Sons);
     273      }
     274        }
     275   gmesh->BuildConnectivity();
     276   
     277}
     278/*}}}*/
     279void AdaptiveMeshRefinement::GetMesh(int &nvertices,int &nelements,double** px,double** py, int** pelements){/*{{{*/
     280
     281        /* IMPORTANT! pelements are in Matlab indexing
     282           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.
     285                Avoid to call this method before Refinement Process.*/
     286
     287        /*Intermediaries */
     288        TPZGeoMesh* gmesh = this->currentmesh;//itapopo confirmar
     289       
     290        long sid,nodeindex;
     291        int nconformelements,nconformvertices; 
     292        int* newelements                        = NULL;
     293        double* newmeshX                        = NULL;//xNew<double>(ntotalvertices);
     294        double* newmeshY                        = NULL;//xNew<double>(ntotalvertices);
     295        TPZGeoEl* geoel                 = NULL;
     296        long* vertex_index2sid  = xNew<long>(gmesh->NNodes());
     297        this->sid2index.clear();
     298       
     299        /*Get mesh coords */
     300        //for(int i=0;i<ntotalvertices;i++ ){
     301        //      TPZVec<REAL> coords(3,0.);
     302        //      gmesh->NodeVec()[i].GetCoordinates(coords);
     303        //      newmeshX[i] = coords[0];
     304        //      newmeshY[i] = coords[1];
     305        //}
     306
     307        /*Fill in the vertex_index2sid vector with non usual index value*/
     308        for(int i=0;i<gmesh->NNodes();i++) vertex_index2sid[i]=-1;
     309       
     310        /*Get elements without sons and fill in the vertex_index2sid with used vertices (indexes) */
     311        sid=0;
     312        for(int i=0;i<gmesh->NElements();i++){//over gmesh elements index
     313                geoel=gmesh->ElementVec()[i];
     314                if(!geoel) continue;
     315                if(geoel->HasSubElement()) continue;
     316                if(geoel->MaterialId() != this->GetElemMaterialID()) continue;
     317                this->sid2index.push_back(i);//keep the element index
     318                for(int j=0;j<this->elementswidth;j++){
     319        nodeindex=geoel->NodeIndex(j);
     320        if(vertex_index2sid[nodeindex]==-1){
     321                vertex_index2sid[nodeindex]=sid;
     322                                sid++;
     323                        }
     324      }
     325        }
     326
     327        nconformelements        = (int)this->sid2index.size();
     328        nconformvertices        = (int)sid;
     329        newelements                     = xNew<int>(nconformelements*this->elementswidth);
     330        newmeshX                                = xNew<double>(nconformvertices);
     331   newmeshY                             = xNew<double>(nconformvertices);
     332
     333        for(int i=0;i<nconformvertices;i++){//over the TPZNode index (fill in the ISSM vertices coords)
     334                sid = vertex_index2sid[i];
     335                if(sid!=-1){
     336                        TPZVec<REAL> coords(3,0.);
     337                        gmesh->NodeVec()[i].GetCoordinates(coords);
     338                        newmeshX[sid] = coords[0];
     339                        newmeshY[sid] = coords[1];
     340                }
     341        }
     342               
     343        for(int i=0;i<this->sid2index.size();i++){//over the sid (fill the ISSM elements)
     344                for(int j=0;j<this->elementswidth;j++) {
     345                        geoel   = gmesh->ElementVec()[this->sid2index[i]];
     346                        sid     = vertex_index2sid[geoel->NodeIndex(j)];
     347                        newelements[i*this->elementswidth+j]=(int)sid+1;//C to Matlab indexing
     348                }
     349        }
     350 
     351        /*Setting outputs*/
     352        nvertices       = nconformvertices;
     353        nelements       = nconformelements;
     354        *px                     = newmeshX;
     355        *py                = newmeshY;
     356        *pelements      = newelements;
     357   
     358        /*Cleanup*/
     359        xDelete<long>(vertex_index2sid);
     360
     361}
     362/*}}}*/
     363void AdaptiveMeshRefinement::FindElements(int &numberofpoints,double* xp,double* yp,TPZGeoMesh *gmesh,int &hlevel,std::vector<int> &elements){/*{{{*/
     364
     365        if(!gmesh) _error_("Impossible to set elements: gmesh is NULL!\n");
     366
     367        if(false){
     368                this->AllElements(gmesh,elements); //uniform, refine all elements!
     369                return;
     370        }
     371
     372        /*Intermediaries*/
     373        elements.clear();
     374        double D1               = this->regionlevel1;
     375        double Dhmax    = this->regionlevelmax;
     376        int hmax                        = this->levelmax;
     377        double alpha    = (hmax==1) ? 0. : log(D1/Dhmax)/(hmax-1.);
     378        double Di               = D1/exp(alpha*(hlevel-1));
     379        int side2D              = 6;
     380        double distance,value;
     381   
     382        /*Find elements near the points */
     383        for(int i=0;i<gmesh->NElements();i++){
     384                if(gmesh->Element(i)->MaterialId()!=this->GetElemMaterialID()) continue;
     385                if(gmesh->Element(i)->HasSubElement()) continue;
     386                if(gmesh->Element(i)->Level()>=hlevel) continue;
     387                TPZVec<REAL> qsi(2,0.);
     388                TPZVec<REAL> centerPoint(3,0.);
     389                gmesh->Element(i)->CenterPoint(side2D, qsi);
     390                gmesh->Element(i)->X(qsi, centerPoint);
     391                distance = Di;
     392                for (int j=0;j<numberofpoints;j++){
     393                        value = std::sqrt( (xp[j]-centerPoint[0])*(xp[j]-centerPoint[0])+(yp[j]-centerPoint[1] )*(yp[j]-centerPoint[1]) );//sqrt( (x2-x1)^2 + (y2-y1)^2 )
     394                        if(value<distance) distance=value; //min distance to the point
     395                } 
     396                if(distance<Di) elements.push_back(i);
     397        }
     398
     399}
     400/*}}}*/
     401void AdaptiveMeshRefinement::AllElements(TPZGeoMesh *gmesh,std::vector<int> &elements){/*{{{*/
    435402    /* Uniform refinement. This refines the entire mesh */
    436403    int nelements = gmesh->NElements();
     404         elements.clear();
    437405    for(int i=0;i<nelements;i++){
    438406        if(gmesh->Element(i)->MaterialId()!=this->GetElemMaterialID()) continue;
    439407        if(gmesh->Element(i)->HasSubElement()) continue;
    440         ElemVec.push_back(i);
     408        elements.push_back(i);
    441409    }
    442410}
    443411/*}}}*/
    444 void AdaptiveMeshRefinement::TagElementsNearGroundingLine(TPZGeoMesh *gmesh,std::vector<TPZVec<REAL> > &GLvec,int &hlevel,std::vector<int> &ElemVec){/*{{{*/
    445    
    446     /* Tag elements near grounding line */
    447          double D1              = this->regionlevel1;
    448          double Dhmax  = this->regionlevelmax;
    449          int hmax               = this->levelmax;
    450     double alpha        = (hmax==1) ? 0. : log(D1/Dhmax)/(hmax-1.);
    451          double Di              = D1/exp(alpha*(hlevel-1));
    452    
    453     for(int i=0;i<gmesh->NElements();i++){
    454        
    455         if(gmesh->Element(i)->MaterialId()!=this->GetElemMaterialID()) continue;
    456         if(gmesh->Element(i)->HasSubElement()) continue;
    457         if(gmesh->Element(i)->Level()>=hlevel) continue;
    458        
    459         const int side2D = 6;
    460         TPZVec<REAL> qsi(2,0.);
    461         TPZVec<REAL> centerPoint(3,0.);
    462         gmesh->Element(i)->CenterPoint(side2D, qsi);
    463         gmesh->Element(i)->X(qsi, centerPoint);
    464        
    465         REAL distance = Di;
    466        
    467         for (int j = 0; j < GLvec.size(); j++) {
    468            
    469             REAL value = ( GLvec[j][0] - centerPoint[0] ) * ( GLvec[j][0] - centerPoint[0] ); // (x2-x1)^2
    470             value += ( GLvec[j][1] - centerPoint[1] ) * ( GLvec[j][1] - centerPoint[1] );// (y2-y1)^2
    471             value = std::sqrt(value); //Radius
    472            
    473             //finding the min distance to the grounding line
    474             if(value < distance) distance = value;
    475            
    476         }
    477        
    478         if(distance < Di) ElemVec.push_back(i);
    479     }
    480    
    481 }
    482 /*}}}*/
    483 void AdaptiveMeshRefinement::CreateInitialMesh(int &nvertices,int &nelements,int &nsegments,int &width,double* x,double* y,double* z,int* elements,int* segments){/*{{{*/
    484 
    485         /*IMPORTANT! elements come in Matlab indexing*/
    486         /*NEOPZ works only in C indexing*/
    487        
    488         _assert_(nvertices>0);
    489    _assert_(nelements>0);
     412void AdaptiveMeshRefinement::CreateInitialMesh(int &nvertices,int &nelements,int &width,double* x,double* y,int* elements){/*{{{*/
     413
     414        /* IMPORTANT! elements come in Matlab indexing
     415                NEOPZ works only in C indexing*/
     416       
     417        if(nvertices<=0) _error_("Impossible to create initial mesh: nvertices is <= 0!\n");
     418   if(nelements<=0) _error_("Impossible to create initial mesh: nelements is <= 0!\n");
    490419        this->SetElementWidth(width);
    491420
    492421    /*Verify and creating initial mesh*/
    493    if(this->fathermesh) _error_("Initial mesh already exists!");
     422   if(this->fathermesh || this->currentmesh) _error_("Initial mesh already exists!");
    494423   
    495424   this->fathermesh = new TPZGeoMesh();
    496         this->fathermesh->NodeVec().Resize( nvertices );
    497 
    498     /*Set the vertices (geometric nodes in NeoPZ context)*/
     425        this->fathermesh->NodeVec().Resize(nvertices);
     426
     427        /*Set the vertices (geometric nodes in NeoPZ context)*/
    499428        for(int i=0;i<nvertices;i++){ 
    500429      /*x,y,z coords*/
    … …  
    502431      coord[0]= x[i];
    503432      coord[1]= y[i];
    504       coord[2]= z[i];
     433      coord[2]= 0.;
    505434      /*Insert in the mesh*/
    506435      this->fathermesh->NodeVec()[i].SetCoord(coord);
    … …  
    512441   const int mat = this->GetElemMaterialID();
    513442   TPZManVector<long> elem(this->elementswidth,0);
    514    
    515         for(int iel=0;iel<nelements;iel++){
    516 
    517                 for(int jel=0;jel<this->elementswidth;jel++) elem[jel]=elements[iel*this->elementswidth+jel]-1;//Convert Matlab to C indexing
    518 
     443   this->sid2index.clear();
     444
     445        for(int i=0;i<nelements;i++){
     446                for(int j=0;j<this->elementswidth;j++) elem[j]=elements[i*this->elementswidth+j]-1;//Convert Matlab to C indexing
    519447      /*reftype = 0: uniform, fast / reftype = 1: uniform and non-uniform (avoid hanging nodes), it is not too fast */
    520       const int reftype = 0;
     448      const int reftype = 1;
    521449      switch(this->elementswidth){
    522                         case 3: this->fathermesh->CreateGeoElement(ETriangle, elem, mat, index, reftype);       break;
    523          case 4: this->fathermesh->CreateGeoElement(ETetraedro, elem, mat, index, reftype); DebugStop(); break;
    524                         case 6: this->fathermesh->CreateGeoElement(EPrisma, elem, mat, index, reftype); DebugStop(); break;
    525          default:       DebugStop();//itapopo _error_("mesh not supported yet");
     450                        case 3: this->fathermesh->CreateGeoElement(ETriangle,elem,mat,index,reftype);   break;
     451         default:       _error_("mesh not supported yet");
    526452                }
    527        
    528453      /*Define the element ID*/       
    529       this->fathermesh->ElementVec()[index]->SetId(iel);
    530         }
    531    
    532    /*Generate the 1D segments elements (boundary)*/
    533    const int matboundary = this->GetBoundaryMaterialID();
    534    TPZManVector<long> boundary(2,0.);
    535    
    536    for(int iel=nelements;iel<nelements+nsegments;iel++){     
    537                 boundary[0] = segments[(iel-nelements)*2+0]-1;//Convert Matlab to C indexing
    538       boundary[1] = segments[(iel-nelements)*2+1]-1;//Convert Matlab to C indexing
    539       /*reftype = 0: uniform, fast / reftype = 1: uniform and non-uniform (avoid hanging nodes), it is not too fast */
    540       const int reftype = 0;
    541       this->fathermesh->CreateGeoElement(EOned, boundary, matboundary, index, reftype);//cria elemento unidimensional
    542       this->fathermesh->ElementVec()[index]->SetId(iel);
    543         }
    544    
     454      this->fathermesh->ElementVec()[index]->SetId(i);
     455                /*Initialize sid2index*/
     456                this->sid2index.push_back((int)index);
     457        }
    545458   /*Build element and node connectivities*/
    546459   this->fathermesh->BuildConnectivity();
    547    
    548         /*Create previous mesh as a copy of father mesh*/
    549    this->previousmesh = new TPZGeoMesh(*this->fathermesh);
    550 
    551 }
    552 /*}}}*/
    553 #include "pzgeotriangle.h" //itapopo
    554 #include "pzreftriangle.h" //itapopo
    555 using namespace pzgeom;
     460        /*Set current mesh*/
     461        this->currentmesh=new TPZGeoMesh(*this->fathermesh);
     462}
     463/*}}}*/
    556464TPZGeoMesh* AdaptiveMeshRefinement::CreateRefPatternMesh(TPZGeoMesh* gmesh){/*{{{*/
    557465       
    … …  
    646554/*}}}*/
    647555void AdaptiveMeshRefinement::SetElementWidth(int &width){/*{{{*/
    648     this->elementswidth = width;
    649 }
    650 /*}}}*/
    651 void AdaptiveMeshRefinement::CheckMesh(int &nvertices,int &nelements,int &nsegments,int &width,double** px,double** py,double** pz,int** pelements, int** psegments){/*{{{*/
    652 
    653     /*Basic verification*/
    654     if( !(nvertices > 0) || !(nelements > 0) ) DebugStop(); //itapopo verificar se irá usar o _assert_
    655    
    656     if ( !(width == 3) && !(width == 4) && !(width == 6) ) DebugStop(); // itapopo verifcar se irá usar o _assert_
    657    
    658     if( !px || !py || !pz || !pelements ) DebugStop(); // itapopo verifcar se irá usar o _assert_
    659    
    660     /*Verify if there are orphan nodes*/
    661     std::set<int> elemvertices;
    662     elemvertices.clear();
    663     for(int i = 0; i < nelements; i++){
    664         for(int j = 0; j < width; j++) {
    665             elemvertices.insert((*pelements)[i*width+j]);
    666                   }
    667          }
    668    
    669     if( elemvertices.size() != nvertices ) DebugStop();//itapopo verificar se irá usar o _assert_
    670        
    671     //Verify if there are inf or NaN in coords
    672     for(int i = 0; i < nvertices; i++) if(isnan((*px)[i]) || isinf((*px)[i])) DebugStop();
    673     for(int i = 0; i < nvertices; i++) if(isnan((*py)[i]) || isinf((*py)[i])) DebugStop();
    674     for(int i = 0; i < nvertices; i++) if(isnan((*pz)[i]) || isinf((*pz)[i])) DebugStop();
    675    
    676          for(int i = 0; i < nelements; i++){
    677         for(int j = 0; j < width; j++){
    678             if( isnan((*pelements)[i*width+j]) || isinf((*pelements)[i*width+j]) ) DebugStop();
    679         }
    680     }
    681    
    682 }
    683 /*}}}*/
     556        if(width!=3) _error_("elementswidth not supported yet!");
     557   this->elementswidth = width;
     558}
     559/*}}}*/
     560void AdaptiveMeshRefinement::CheckMesh(int &nvertices,int &nelements,int &width,double** px,double** py,int** pelements){/*{{{*/
     561
     562        /*Basic verification*/
     563        if(nvertices<=0) _error_("Impossible to continue: nvertices <=0!\n");
     564        if(nelements<=0) _error_("Impossible to continue: nelements <=0!\n");
     565        if(width!=3) _error_("Impossible to continue: width !=3!\n");
     566        if(!px) _error_("Impossible to continue: px is NULL!\n");
     567        if(!py) _error_("Impossible to continue: py is NULL!\n");
     568        if(!pelements) _error_("Impossible to continue: pelements is NULL!\n");
     569
     570        /*Verify if there are orphan nodes*/
     571        std::set<int> elemvertices;
     572        elemvertices.clear();
     573        for(int i=0;i<nelements;i++){
     574                for(int j=0;j<width;j++) {
     575                        elemvertices.insert((*pelements)[i*width+j]);
     576                }
     577        }
     578        if(elemvertices.size()!=nvertices) _error_("Impossible to continue: elemvertices.size() != nvertices!\n");
     579       
     580        //Verify if there are inf or NaN in coords
     581        for(int i=0;i<nvertices;i++){
     582                if(isnan((*px)[i]) || isinf((*px)[i])) _error_("Impossible to continue: px i=" << i <<" is NaN or Inf!\n");
     583                if(isnan((*py)[i]) || isinf((*py)[i])) _error_("Impossible to continue: py i=" << i <<" is NaN or Inf!\n");
     584        }
     585        for(int i=0;i<nelements;i++){
     586                for(int j=0;j<width;j++){
     587                        if(isnan((*pelements)[i*width+j]) || isinf((*pelements)[i*width+j]) ) _error_("Impossible to continue: px i=" << i <<" is NaN or Inf!\n");
     588                }
     589        }
     590   
     591}
     592/*}}}*/

  • issm/trunk-jpl/src/c/classes/AdaptiveMeshRefinement.h

    r21672 r21806  
    4545        /*Public methods*/
    4646        /* Constructor, destructor etc*/
    47         AdaptiveMeshRefinement();                                                                                                                               // Default constructor
    48         AdaptiveMeshRefinement(const AdaptiveMeshRefinement &cp);                                       // Copy constructor
    49         AdaptiveMeshRefinement & operator= (const AdaptiveMeshRefinement &cp);  // Operator of copy
    50         virtual ~AdaptiveMeshRefinement();                                                                                                      // Destructor
     47        AdaptiveMeshRefinement();                                                                                                                       
     48        AdaptiveMeshRefinement(const AdaptiveMeshRefinement &cp);                                       
     49        AdaptiveMeshRefinement & operator= (const AdaptiveMeshRefinement &cp); 
     50        virtual ~AdaptiveMeshRefinement();                                                                                             
    5151
    5252    /*Savable methods*/
    53         virtual int ClassId() const;                                            // ClassId to save the class
    54    virtual void Read(TPZStream &buf, void *context);                                                            // Read this class
    55    virtual void Write(TPZStream &buf, int withclassid);                    // Write this class, using ClassId to identify
     53        virtual int ClassId() const;                                   
     54   virtual void Read(TPZStream &buf,void *context);                                                     
     55   virtual void Write(TPZStream &buf,int withclassid);
    5656   
    5757        /*General methods*/
    58         void CleanUp();                                                                                                                                                 // Clean all attributes
    59         void Initialize();                                                                                                                                              // Initialize the attributes with NULL and values out of usually range
    60         void SetLevelMax(int &h);                                               // Define the max level of refinement
    61    void SetRegions(double &D1,double Dhmax);                                                                            // Define the regions which will be refined
    62         void SetElementWidth(int &width);                                       // Define elements width
    63         void ExecuteRefinement(int &type_process,double *vx,double *vy,double *masklevelset,int &nvertices,int &nelements,int &nsegments,double** px,double** py,double** pz,int** pelements,int** psegments=NULL);                                     // A new mesh will be created and refined. This returns the new mesh
    64         void CreateInitialMesh(int &nvertices,int &nelements,int &nsegments,int &width,double* x,double* y,double* z,int* elements,int* segments=NULL); // Create a NeoPZ geometric mesh by coords and elements
     58        void CleanUp();
     59        void Initialize();
     60        void SetLevelMax(int &h);
     61   void SetRegions(double &D1,double Dhmax);
     62        void SetElementWidth(int &width);
     63        void Execute(bool &amr_verbose,int &numberofelements,
     64                                                double* partiallyfloatedelements,double *masklevelset,double* deviatorictensorerror,double* thicknesserror,
     65                                                int &newnumberofvertices,int &newnumberofelements,double** x,double** y,int** elementslist);
     66        void CreateInitialMesh(int &nvertices,int &nelements,int &width,double* x,double* y,int* elements);
    6567        TPZGeoMesh* CreateRefPatternMesh(TPZGeoMesh* gmesh);
    66         void CheckMesh(int &nvertices,int &nelements,int &nsegments,int &width,double** px,double** py,double** pz,int** pelements,int** psegments=NULL); // Check the consistency of the mesh
     68        void CheckMesh(int &nvertices,int &nelements,int &width,double** px,double** py,int** pelements);
    6769
    6870private:
    69 
    7071        /*Private attributes*/
    71    int elementswidth;                                                      // Geometric nodes for element: 3 == Tria, 4 == Tetra, 6 == Penta
    72    int levelmax;                                                           // Max level of refinement
    73         double regionlevel1;                                                                                                                                            // Region which will be refined with level 1
    74         double regionlevelmax;                                                                                                                                  // Region which will be refined with level max
    75         TPZGeoMesh *fathermesh;                                                                                                                                 // Father Mesh is the entire mesh without refinement
    76         TPZGeoMesh *previousmesh;                                                                                                                               // Previous mesh is a refined mesh of last step
     72   int elementswidth;                                    // Geometric nodes for element: 3 == Tria, 4 == Tetra, 6 == Penta
     73   int levelmax;                                         // Max level of refinement
     74        double regionlevel1;                                                                                            // Region which will be refined with level 1
     75        double regionlevelmax;                                                                                  // Region which will be refined with level max
     76        std::vector<int> sid2index;                                                                     // Vector that keeps index of PZGeoMesh elements used in the ISSM mesh (sid)
     77        TPZGeoMesh *fathermesh;                                                                                 // Father Mesh is the entire mesh without refinement
     78        TPZGeoMesh *currentmesh;                                                                                // Current Mesh is the refined mesh
    7779
    7880        /*Private methods*/
    79    void RefinementProcess(TPZGeoMesh *gmesh,std::vector<TPZVec<REAL> > &GLvec);  // Start the refinement process
    80         void RefineMesh(TPZGeoMesh *gmesh, std::vector<int> &ElemVec);                                  // Refine the elements in ElemVec
    81    void RefineMeshToAvoidHangingNodes(TPZGeoMesh *gmesh);                        // Refine the elements to avoid hanging nodes
    82         void SetElementsToRefine(TPZGeoMesh *gmesh,std::vector<TPZVec<REAL> > &GLvec,int &hlevel, std::vector<int> &ElemVec);   //Define wich elements will be refined
    83    void TagAllElements(TPZGeoMesh *gmesh,std::vector<int> &ElemVec);                             // This tag all elements to be refined, that is, refine all elements
    84    void TagElementsNearGroundingLine(TPZGeoMesh *gmesh,std::vector<TPZVec<REAL> > &GLvec,int &hlevel,std::vector<int> &ElemVec);    // This tag elements near the grounding line
    85    void CalcGroundingLinePosition(double *masklevelset,std::vector<TPZVec<REAL> > &GLvec);      // Calculate the grounding line position using previous mesh
    86         void GetMesh(TPZGeoMesh *gmesh,int &nvertices,int &nelements,int &nsegments,double** px,double** py,double** pz,int** pelements,int** psegments=NULL); // Return coords and elements in ISSM data structure
    87    inline int GetElemMaterialID(){return 1;}                               // Return element material ID
    88    inline int GetBoundaryMaterialID(){return 2;}                           // Return segment (2D boundary) material ID
     81   void RefinementProcess(bool &amr_verbose,double* partiallyfloatedelements,double* masklevelset,double* deviatorictensorerror,double* thicknesserror);
     82        void RefineMesh(TPZGeoMesh *gmesh,std::vector<int> &elements);
     83   void RefineMeshToAvoidHangingNodes(TPZGeoMesh *gmesh);
     84        void GetMesh(int &nvertices,int &nelements,double** px,double** py,int** pelements);
     85        void FindElements(int &numberofpoints,double* xp,double* yp,TPZGeoMesh *gmesh,int &hlevel,std::vector<int> &elements);
     86   void AllElements(TPZGeoMesh *gmesh,std::vector<int> &elements);
     87   inline int GetElemMaterialID(){return 1;}         
    8988};
    9089

  • issm/trunk-jpl/src/c/classes/FemModel.cpp

    r21803 r21806  
    26932693void FemModel::WriteMeshInResults(void){/*{{{*/
    26942694
     2695        //itapopo
     2696        #ifdef _HAVE_NEOPZ_
     2697        this->WriteErrorEstimatorsInResults();
     2698        #endif
     2699        //itapopo
     2700
    26952701        int step                                        = -1;
    26962702        int numberofelements = -1;
    … …  
    27512757        xDelete<IssmDouble>(z);
    27522758        xDelete<int>(elementslist);
     2759}
     2760/*}}}*/
     2761void FemModel::WriteErrorEstimatorsInResults(void){/*{{{*/
     2762
     2763   int step                   = -1;
     2764   int numberofelements       = -1;
     2765   IssmDouble time            = -1;
     2766   IssmDouble* stresserror    = NULL;
     2767   IssmDouble* thicknesserror = NULL;
     2768
     2769   if(!this->elements || !this->vertices || !this->results || !this->parameters) return;
     2770
     2771   parameters->FindParam(&step,StepEnum);
     2772   parameters->FindParam(&time,TimeEnum);
     2773   numberofelements=this->elements->NumberOfElements();
     2774
     2775   /*Compute the deviatoric stress tensor*/
     2776   this->ZZErrorEstimator(&stresserror);
     2777
     2778   /*Compute the thickness error*/
     2779   this->ThicknessZZErrorEstimator(&thicknesserror);
     2780
     2781   /*Write error estimators in Results*/
     2782   this->results->AddResult(new GenericExternalResult<IssmDouble*>(this->results->Size()+1,DeviatoricStressErrorEstimatorEnum,
     2783                                                                  stresserror,numberofelements,1,step,time));
     2784
     2785   this->results->AddResult(new GenericExternalResult<IssmDouble*>(this->results->Size()+1,ThicknessErrorEstimatorEnum,
     2786                                                                  thicknesserror,numberofelements,1,step,time));
     2787   /*Cleanup*/
     2788   xDelete<IssmDouble>(stresserror);
     2789   xDelete<IssmDouble>(thicknesserror);
     2790
     2791   return;
    27532792}
    27542793/*}}}*/
    … …  
    32123251void FemModel::SmoothedDeviatoricStressTensor(IssmDouble** ptauxx,IssmDouble** ptauyy,IssmDouble** ptauxy){/*{{{*/
    32133252       
    3214         this->DeviatoricStressx();//itapopo
    3215 
    32163253        int elementswidth                                                       = this->GetElementsWidth();//just 2D mesh, tria elements
    32173254   int numberofvertices                                         = this->vertices->NumberOfVertices();
    3218 
    32193255   IssmDouble weight                                            = 0.;
    32203256        IssmDouble*     tauxx                                                   = NULL;
    … …  
    32313267   Vector<IssmDouble>* vectotalweight   = new Vector<IssmDouble>(numberofvertices);
    32323268       
    3233    for(int i=0;i<this->elements->Size();i++){
     3269        /*Update the Deviatoric Stress tensor over the elements*/
     3270        this->DeviatoricStressx();
     3271       
     3272   /*Calculate the Smoothed Deviatoric Stress tensor*/
     3273        for(int i=0;i<this->elements->Size();i++){
    32343274      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
    32353275      element->GetInputListOnVertices(deviatoricstressxx,DeviatoricStressxxEnum);
    … …  
    32993339void FemModel::ZZErrorEstimator(IssmDouble** pelementerror){/*{{{*/
    33003340
    3301         /*Compute the Zienkiewicz and Zhu (ZZ) error estimator.
     3341        /*Compute the Zienkiewicz and Zhu (ZZ) error estimator for the deviatoric stress tensor.
    33023342         * Ref.: Zienkiewicz and Zhu, A Simple Error Estimator and Adaptive Procedure for Practical Engineering Analysis, Int. J. Numer. Meth. Eng, 1987*/
    33033343
    … …  
    33423382                                ftxy+=(tauxy[n]-smoothedtauxy[elem_vertices[n]])*basis[n];
    33433383                        }
    3344                         error+=Jdet*gauss->weight*( std::pow(ftxx,2)+std::pow(ftyy,2)+std::pow(ftxy,2) );
     3384                        error+=Jdet*gauss->weight*( std::pow(ftxx,2)+std::pow(ftyy,2)+std::pow(ftxy,2) ); //e^2
    33453385                }
    33463386                /*Set the error in the global vector*/ 
    33473387      sid=element->Sid();
    3348                 velementerror->SetValue(sid,error,INS_VAL);     
     3388                velementerror->SetValue(sid,std::sqrt(error),INS_VAL);//sqrt( e^2 )
    33493389                /*Cleanup intermediaries*/
    33503390                xDelete<IssmDouble>(xyz_list);
    … …  
    33683408        xDelete<int>(elem_vertices);
    33693409        delete velementerror;
     3410}
     3411/*}}}*/
     3412void FemModel::SmoothedGradThickness(IssmDouble** pdHdx,IssmDouble** pdHdy){/*{{{*/
     3413
     3414   int elementswidth                   = this->GetElementsWidth();//just 2D mesh, tria elements
     3415   int numberofvertices                = this->vertices->NumberOfVertices();
     3416
     3417   IssmDouble weight                   = 0.;
     3418   IssmDouble* dHdx                    = NULL;
     3419   IssmDouble* dHdy                    = NULL;
     3420   IssmDouble* totalweight             = NULL;
     3421   IssmDouble* xyz_list                = NULL;
     3422   IssmDouble* H                       = xNew<IssmDouble>(elementswidth);
     3423   IssmDouble* GradH                   = xNew<IssmDouble>(2);
     3424   int* elem_vertices                  = xNew<int>(elementswidth);
     3425   Vector<IssmDouble>* vecdHdx         = new Vector<IssmDouble>(numberofvertices);
     3426   Vector<IssmDouble>* vecdHdy         = new Vector<IssmDouble>(numberofvertices);
     3427   Vector<IssmDouble>* vectotalweight  = new Vector<IssmDouble>(numberofvertices);
     3428
     3429   for(int i=0;i<this->elements->Size();i++){
     3430      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
     3431      element->GetInputListOnVertices(H,ThicknessEnum);
     3432      element->GetVerticesSidList(elem_vertices);
     3433      element->GetVerticesCoordinates(&xyz_list);
     3434
     3435      /*Get the gradient of thickness at the center point (in fact, GradH is constante over the element)*/
     3436      Gauss* gauss=element->NewGauss(1);
     3437      gauss->GaussPoint(gauss->begin());
     3438      element->ValueP1DerivativesOnGauss(GradH,H,xyz_list,gauss);
     3439
     3440      /*weight to calculate the smoothed grad H*/
     3441      Tria* triaelement = xDynamicCast<Tria*>(element);
     3442      weight            = triaelement->GetArea();//the tria area is a choice for the weight
     3443     
     3444                /*dH/dx*/
     3445      vecdHdx->SetValue(elem_vertices[0],weight*GradH[0],ADD_VAL);
     3446      vecdHdx->SetValue(elem_vertices[1],weight*GradH[0],ADD_VAL);
     3447      vecdHdx->SetValue(elem_vertices[2],weight*GradH[0],ADD_VAL);
     3448      /*dH/dy*/
     3449      vecdHdy->SetValue(elem_vertices[0],weight*GradH[1],ADD_VAL);
     3450      vecdHdy->SetValue(elem_vertices[1],weight*GradH[1],ADD_VAL);
     3451      vecdHdy->SetValue(elem_vertices[2],weight*GradH[1],ADD_VAL);
     3452      /*total weight*/
     3453      vectotalweight->SetValue(elem_vertices[0],weight,ADD_VAL);
     3454      vectotalweight->SetValue(elem_vertices[1],weight,ADD_VAL);
     3455      vectotalweight->SetValue(elem_vertices[2],weight,ADD_VAL);
     3456      /*Cleanup intermediaries*/
     3457      xDelete<IssmDouble>(xyz_list);
     3458      delete gauss;
     3459   }
     3460
     3461   /*Assemble*/
     3462   vecdHdx->Assemble();
     3463   vecdHdy->Assemble();
     3464   vectotalweight->Assemble();
     3465
     3466   /*Serialize*/
     3467   dHdx        = vecdHdx->ToMPISerial();
     3468   dHdy        = vecdHdy->ToMPISerial();
     3469   totalweight = vectotalweight->ToMPISerial();
     3470
     3471   /*Divide for the total weight*/
     3472   for(int i=0;i<numberofvertices;i++){
     3473      _assert_(totalweight[i]>0);
     3474      dHdx[i] = dHdx[i]/totalweight[i];
     3475      dHdy[i] = dHdy[i]/totalweight[i];
     3476   }
     3477
     3478   /*Set output*/
     3479   (*pdHdx) = dHdx;
     3480   (*pdHdy) = dHdy;
     3481
     3482        /*Cleanup*/
     3483   delete vecdHdx;
     3484   delete vecdHdy;
     3485   delete vectotalweight;
     3486   xDelete<IssmDouble>(H);
     3487   xDelete<IssmDouble>(GradH);
     3488   xDelete<IssmDouble>(totalweight);
     3489   xDelete<int>(elem_vertices);
     3490}
     3491/*}}}*/
     3492void FemModel::ThicknessZZErrorEstimator(IssmDouble** pelementerror){/*{{{*/
     3493   /*Compute the Zienkiewicz and Zhu (ZZ) error estimator for the thickness
     3494    * Ref.: Zienkiewicz and Zhu, A Simple Error Estimator and Adaptive Procedure for Practical Engineering Analysis, Int. J. Numer. Meth. Eng, 1987*/
     3495
     3496   IssmDouble Jdet,error,fdHdx,fdHdy;
     3497   int sid;
     3498   int numnodes                     = this->GetElementsWidth();//just 2D mesh, tria elements, P1
     3499   int numberofelements             = this->elements->NumberOfElements();
     3500   IssmDouble* xyz_list             = NULL;
     3501   IssmDouble* smoothed_dHdx        = NULL;
     3502   IssmDouble* smoothed_dHdy        = NULL;
     3503   IssmDouble* H                    = xNew<IssmDouble>(numnodes);
     3504   IssmDouble* GradH                = xNew<IssmDouble>(2);
     3505   IssmDouble* basis                = xNew<IssmDouble>(numnodes);
     3506   int* elem_vertices               = xNew<int>(numnodes);
     3507   Vector<IssmDouble>* velementerror= new Vector<IssmDouble>(numberofelements);
     3508
     3509   /*Get smoothed deviatoric stress tensor*/
     3510   this->SmoothedGradThickness(&smoothed_dHdx,&smoothed_dHdy);
     3511   
     3512        /*Integrate the error over elements*/
     3513   for(int i=0;i<this->elements->Size();i++){
     3514      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
     3515      element->GetInputListOnVertices(H,ThicknessEnum);
     3516      element->GetVerticesSidList(elem_vertices);
     3517      element->GetVerticesCoordinates(&xyz_list);
     3518      /*Get the gradient of thickness*/
     3519      Gauss* gaussH=element->NewGauss(1);
     3520      gaussH->GaussPoint(gaussH->begin());
     3521      element->ValueP1DerivativesOnGauss(GradH,H,xyz_list,gaussH);
     3522      /*Integrate*/
     3523      Gauss* gauss=element->NewGauss(2);
     3524      error=0.;
     3525      for(int ig=gauss->begin();ig<gauss->end();ig++){
     3526         gauss->GaussPoint(ig);
     3527         element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     3528         element->NodalFunctions(basis,gauss);
     3529         fdHdx=0;fdHdy=0;
     3530         for(int n=0;n<numnodes;n++) {
     3531            fdHdx+=(GradH[0]-smoothed_dHdx[elem_vertices[n]])*basis[n];
     3532            fdHdy+=(GradH[1]-smoothed_dHdy[elem_vertices[n]])*basis[n];
     3533         }
     3534         error+=Jdet*gauss->weight*( std::pow(fdHdx,2)+std::pow(fdHdy,2) ); //e^2
     3535      }
     3536      /*Set the error in the global vector*/
     3537      sid=element->Sid();
     3538      velementerror->SetValue(sid,std::sqrt(error),INS_VAL);//sqrt( e^2 )
     3539      /*Cleanup intermediaries*/
     3540      xDelete<IssmDouble>(xyz_list);
     3541      delete gaussH;
     3542      delete gauss;
     3543   }
     3544
     3545   /*Assemble*/
     3546   velementerror->Assemble();
     3547
     3548   /*Serialize and set output*/
     3549   (*pelementerror)=velementerror->ToMPISerial();
     3550
     3551   /*Cleanup*/
     3552   xDelete<IssmDouble>(smoothed_dHdx);
     3553   xDelete<IssmDouble>(smoothed_dHdy);
     3554   xDelete<IssmDouble>(H);
     3555   xDelete<IssmDouble>(GradH);
     3556   xDelete<IssmDouble>(basis);
     3557   xDelete<int>(elem_vertices);
     3558   delete velementerror;
     3559}
     3560/*}}}*/
     3561void FemModel::MeanGroundedIceLevelSet(IssmDouble** pmasklevelset){/*{{{*/
     3562
     3563   int elementswidth                   = this->GetElementsWidth();
     3564   int numberofelements                = this->elements->NumberOfElements();
     3565   IssmDouble* elementlevelset         = xNew<IssmDouble>(elementswidth);
     3566   Vector<IssmDouble>* vmasklevelset   = new Vector<IssmDouble>(numberofelements);
     3567
     3568   for(int i=0;i<this->elements->Size();i++){
     3569      Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
     3570      element->GetInputListOnVertices(elementlevelset,MaskGroundediceLevelsetEnum);
     3571      int sid = element->Sid();
     3572      vmasklevelset->SetValue(sid,(elementlevelset[0]+elementlevelset[1]+elementlevelset[2])/3.,INS_VAL);
     3573   }
     3574
     3575   /*Assemble*/
     3576   vmasklevelset->Assemble();
     3577
     3578   /*Serialize and set output*/
     3579   (*pmasklevelset)=vmasklevelset->ToMPISerial();
     3580
     3581   /*Cleanup*/
     3582   xDelete<IssmDouble>(elementlevelset);
     3583   delete vmasklevelset;
     3584
    33703585}
    33713586/*}}}*/
    … …  
    41604375
    41614376#ifdef _HAVE_NEOPZ_
    4162 void FemModel::ReMeshNeopz(int* pnumberofvertices,int* pnumberofelements,IssmDouble** px,IssmDouble** py,IssmDouble** pz,int** pelementslist){/*{{{*/
    4163 
    4164         /*elements is in Matlab indexing*/
    4165         int my_rank                                             = IssmComm::GetRank();
    4166         int numberofsegments                    = -1;
    4167         int numberofvertices,numberofelements;
    4168         IssmDouble* vx                                  = NULL; //itapopo this is not being used
    4169         IssmDouble* vy                                  = NULL; //itapopo this is not being used
    4170         IssmDouble* x                                   = NULL;
    4171         IssmDouble* y                                   = NULL;
    4172         IssmDouble* z                                   = NULL;
    4173         int* elementslist                               = NULL;
    4174         int* segments                                   = NULL;
    4175         IssmDouble* masklevelset        = NULL;
    4176         IssmDouble* pelementerror       = NULL;
    4177         const int elementswidth         = this->GetElementsWidth();//just 2D mesh, tria elements
    4178 
    4179         /*Solutions which will be used to refine the elements*/
    4180         this->GetGroundediceLevelSet(&masklevelset);//itapopo verificar se já existe um método igual a esse
    4181 
    4182         //Compute the ZZ error estimator per element
    4183         this->ZZErrorEstimator(&pelementerror);
    4184 
    4185         _printf0_("P Element error\n");
    4186         for(int i=0;i<this->elements->NumberOfElements();i++)   _printf0_(""<<pelementerror[i]<< "\n");
    4187         _printf0_("\n");
     4377void FemModel::ReMeshNeopz(int* pnewnumberofvertices,int* pnewnumberofelements,IssmDouble** pnewx,IssmDouble** pnewy,IssmDouble** pnewz,int** pnewelementslist){/*{{{*/
     4378
     4379        /*pnewelementslist keep vertices in Matlab indexing*/
     4380   int my_rank                         = IssmComm::GetRank();
     4381   bool amr_verbose                    = true; //itapopo
     4382   IssmDouble* x                       = NULL;
     4383   IssmDouble* y                       = NULL;
     4384   IssmDouble* z                       = NULL;
     4385   int* elementslist                   = NULL;
     4386   int oldnumberofelements             = this->elements->NumberOfElements();
     4387   int newnumberofvertices                              = -1;
     4388        int newnumberofelements                                 = -1;
     4389        IssmDouble* partiallyfloatedelements= NULL;//itapopo verify if it will be used
     4390   IssmDouble* masklevelset            = NULL;
     4391   IssmDouble* deviatorictensorerror   = NULL;
     4392   IssmDouble* thicknesserror          = NULL; 
     4393
     4394        /*Get the elements in which grounding line goes through*/
     4395   //itapopo verificar se irá usar esse this->GetPartiallyFloatedElements(&partiallyfloatedelements);
     4396   /*Get mean mask level set over each element*/
     4397   this->MeanGroundedIceLevelSet(&masklevelset);
     4398   /*Get the deviatoric error estimator*/
     4399   this->ZZErrorEstimator(&deviatorictensorerror);
     4400   /*Get the thickness error estimator*/
     4401   this->ThicknessZZErrorEstimator(&thicknesserror);
    41884402
    41894403        if(my_rank==0){
    4190                 int type_process=1; //1: it refines father mesh. See AdaptiveMeshRefinement.h (.cpp)
    4191                 this->amr->ExecuteRefinement(type_process,vx,vy,masklevelset,
    4192                                         numberofvertices,numberofelements,numberofsegments,&x,&y,&z,&elementslist,&segments);
    4193                 if(numberofvertices<=0 || numberofelements<=0 /*|| newnumberofsegments<=0*/) _error_("Error in the refinement process.");
     4404                this->amr->Execute(amr_verbose,oldnumberofelements,partiallyfloatedelements,masklevelset,deviatorictensorerror,thicknesserror,
     4405                           newnumberofvertices,newnumberofelements,&x,&y,&elementslist);
     4406      z=xNewZeroInit<IssmDouble>(newnumberofvertices);
     4407                if(newnumberofvertices<=0 || newnumberofelements<=0) _error_("Error in the ReMeshNeopz.");
    41944408        }
    41954409        else{
    4196                 x=xNew<IssmDouble>(numberofvertices);
    4197                 y=xNew<IssmDouble>(numberofvertices);
    4198                 z=xNew<IssmDouble>(numberofvertices);
    4199                 elementslist=xNew<int>(numberofelements*this->GetElementsWidth());
     4410                x=xNew<IssmDouble>(newnumberofvertices);
     4411                y=xNew<IssmDouble>(newnumberofvertices);
     4412                z=xNew<IssmDouble>(newnumberofvertices);
     4413                elementslist=xNew<int>(newnumberofelements*this->GetElementsWidth());
    42004414        }
    42014415
    42024416        /*Send new mesh to others CPU*/
    4203         ISSM_MPI_Bcast(&numberofvertices,1,ISSM_MPI_INT,0,IssmComm::GetComm());
    4204         ISSM_MPI_Bcast(&numberofelements,1,ISSM_MPI_INT,0,IssmComm::GetComm());
    4205         ISSM_MPI_Bcast(x,numberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());       
    4206         ISSM_MPI_Bcast(y,numberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());       
    4207         ISSM_MPI_Bcast(z,numberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());       
    4208         ISSM_MPI_Bcast(elementslist,numberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm());     
     4417        ISSM_MPI_Bcast(&newnumberofvertices,1,ISSM_MPI_INT,0,IssmComm::GetComm());
     4418        ISSM_MPI_Bcast(&newnumberofelements,1,ISSM_MPI_INT,0,IssmComm::GetComm());
     4419        ISSM_MPI_Bcast(x,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
     4420        ISSM_MPI_Bcast(y,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
     4421        ISSM_MPI_Bcast(z,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
     4422        ISSM_MPI_Bcast(elementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm());   
    42094423
    42104424        /*Assign the pointers*/
    4211         (*pelementslist) = elementslist; //Matlab indexing
    4212         (*px)                             = x;
    4213         (*py)                             = y;
    4214         (*pz)                             = z;
    4215         *pnumberofelements = numberofelements;
    4216         *pnumberofvertices = numberofvertices;
     4425        (*pnewelementslist)     = elementslist; //Matlab indexing
     4426        (*pnewx)                                        = x;
     4427        (*pnewy)                                        = y;
     4428        (*pnewz)                                        = z;
     4429        *pnewnumberofvertices= newnumberofvertices;
     4430        *pnewnumberofelements= newnumberofelements;
    42174431
    42184432        /*Cleanup*/
    4219         if(segments) xDelete<int>(segments);
    42204433        xDelete<IssmDouble>(masklevelset);
    4221         xDelete<IssmDouble>(pelementerror);
     4434        xDelete<IssmDouble>(deviatorictensorerror);
     4435   xDelete<IssmDouble>(thicknesserror);
    42224436
    42234437}
    … …  
    42294443        int numberofvertices                    = this->vertices->NumberOfVertices();
    42304444        int numberofelements                    = this->elements->NumberOfElements();
    4231         int numberofsegments                    = 0; //used on matlab
    42324445        IssmDouble* x                                   = NULL;
    42334446        IssmDouble* y                                   = NULL;
    … …  
    42824495                        //this->amr->SetLevelMax(levelmax); //Set max level of refinement
    42834496                        //this->amr->SetRegions(regionlevel1,regionlevelmax);
    4284                         this->amr->CreateInitialMesh(numberofvertices,numberofelements,numberofsegments,elementswidth,x,y,z,elements,NULL);
     4497                        this->amr->CreateInitialMesh(numberofvertices,numberofelements,elementswidth,x,y,elements);
    42854498                }
    42864499                this->amr->SetLevelMax(levelmax); //Set max level of refinement
    … …  
    42984511
    42994512   /*Initialize the global variable of refinement patterns*/
    4300    gRefDBase.InitializeUniformRefPattern(EOned);
    43014513   gRefDBase.InitializeUniformRefPattern(ETriangle);
    43024514
    4303     //gRefDBase.InitializeRefPatterns();
    43044515   /*Insert specifics patterns to ISSM core*/
    43054516   std::string filepath  = REFPATTERNDIR;

  • issm/trunk-jpl/src/c/classes/FemModel.h

    r21802 r21806  
    170170                void ElementsAndVerticesPartitioning(int& newnumberofvertices,int& newnumberofelements,int& elementswidth,int* newelementslist,bool** pmy_elements,int** pmy_vertices);
    171171                void WriteMeshInResults(void);
    172                 void SetRefPatterns(void);
     172                void WriteErrorEstimatorsInResults(void);
    173173                void SmoothedDeviatoricStressTensor(IssmDouble** ptauxx,IssmDouble** ptauyy,IssmDouble** ptauxy); //nodal values, just for SSA-P1: TauXX, TauYY, TauXY
    174174                void ZZErrorEstimator(IssmDouble** pelementerror);
     175                void SmoothedGradThickness(IssmDouble** pdHdx,IssmDouble** pdHdy);
     176                void ThicknessZZErrorEstimator(IssmDouble** pelementerror);
     177                void MeanGroundedIceLevelSet(IssmDouble** pmasklevelset);
    175178
    176179                #ifdef _HAVE_BAMG_
    … …  
    183186                void ReMeshNeopz(int* pnewnumberofvertices,int* pnewnumberofelements,IssmDouble** pnewx,IssmDouble** pnewy,IssmDouble** pnewz,int** pnewelementslist);
    184187                void InitializeAdaptiveRefinementNeopz(void);
     188                void SetRefPatterns(void);
    185189                #endif
    186190};

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r21802 r21806  
    844844        AmrFieldEnum,
    845845        AmrErrEnum,
     846        DeviatoricStressErrorEstimatorEnum,
     847        ThicknessErrorEstimatorEnum,
    846848        /*}}}*/
    847849        ParametersENDEnum,

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r21802 r21806  
    820820                case AmrFieldEnum : return "AmrField";
    821821                case AmrErrEnum : return "AmrErr";
     822                case DeviatoricStressErrorEstimatorEnum : return "DeviatoricStressErrorEstimator";
     823                case ThicknessErrorEstimatorEnum : return "ThicknessErrorEstimator";
    822824                case ParametersENDEnum : return "ParametersEND";
    823825                case XYEnum : return "XY";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r21802 r21806  
    838838              else if (strcmp(name,"AmrField")==0) return AmrFieldEnum;
    839839              else if (strcmp(name,"AmrErr")==0) return AmrErrEnum;
     840              else if (strcmp(name,"DeviatoricStressErrorEstimator")==0) return DeviatoricStressErrorEstimatorEnum;
     841              else if (strcmp(name,"ThicknessErrorEstimator")==0) return ThicknessErrorEstimatorEnum;
    840842              else if (strcmp(name,"ParametersEND")==0) return ParametersENDEnum;
    841843              else if (strcmp(name,"XY")==0) return XYEnum;
    … …  
    873875              else if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum;
    874876              else if (strcmp(name,"Param")==0) return ParamEnum;
    875               else if (strcmp(name,"Moulin")==0) return MoulinEnum;
    876               else if (strcmp(name,"Pengrid")==0) return PengridEnum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"Penpair")==0) return PenpairEnum;
     880              if (strcmp(name,"Moulin")==0) return MoulinEnum;
     881              else if (strcmp(name,"Pengrid")==0) return PengridEnum;
     882              else if (strcmp(name,"Penpair")==0) return PenpairEnum;
    881883              else if (strcmp(name,"Profiler")==0) return ProfilerEnum;
    882884              else if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum;
    … …  
    996998              else if (strcmp(name,"IceVolumeAboveFloatation")==0) return IceVolumeAboveFloatationEnum;
    997999              else if (strcmp(name,"TotalFloatingBmb")==0) return TotalFloatingBmbEnum;
    998               else if (strcmp(name,"TotalGroundedBmb")==0) return TotalGroundedBmbEnum;
    999               else if (strcmp(name,"TotalSmb")==0) return TotalSmbEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"P0")==0) return P0Enum;
     1003              if (strcmp(name,"TotalGroundedBmb")==0) return TotalGroundedBmbEnum;
     1004              else if (strcmp(name,"TotalSmb")==0) return TotalSmbEnum;
     1005              else if (strcmp(name,"P0")==0) return P0Enum;
    10041006              else if (strcmp(name,"P0Array")==0) return P0ArrayEnum;
    10051007              else if (strcmp(name,"P1")==0) return P1Enum;

  • issm/trunk-jpl/src/m/contrib/tsantos/AMRexportVTK.m

    r21678 r21806  
    3737end
    3838
    39 %this is the result structure
     39%this is the result structure (just the Transient solution)
    4040res_struct=model.results;
    4141%checking for results
    … …  
    5151                if(size(sol_struct{i},2)>num_of_timesteps);
    5252                        num_of_timesteps=size(sol_struct{i},2);
    53       outstep=model.timestepping.time_step*model.settings.output_frequency;
     53                        outstep=model.timestepping.time_step*model.settings.output_frequency;
    5454          end
    5555  end
    … …  
    101101                        else
    102102                                timestep = size(sol_struct{j},2);
    103             end
     103                        end
    104104                       
    105105                        %getting the number of fields in the solution
    … …  
    120120                                        s='%e\n';
    121121                                        fprintf(fid,s,sol_struct{j}(timestep).(fieldnames{k}));
    122                     end         
    123             end
     122                                end
     123                        end
     124                        fprintf(fid,'CELL_DATA %s \n',num2str(num_of_elt));
     125                        for k=1:num_of_fields
     126                                if ((numel(sol_struct{j}(timestep).(fieldnames{k})))==num_of_elt);
     127                                        %paraview does not like NaN, replacing
     128                                        nanval=find(isnan(sol_struct{j}(timestep).(fieldnames{k})));
     129                                        sol_struct{j}(timestep).(fieldnames{k})(nanval)=-9999;
     130                                        %also checking for verry small value that mess up
     131                                        smallval=(abs(sol_struct{j}(timestep).(fieldnames{k}))<1.0e-20);
     132                                        sol_struct{j}(timestep).(fieldnames{k})(smallval)=0.0;
     133                                        fprintf(fid,'SCALARS %s float 1 \n',fieldnames{k});
     134                                        fprintf(fid,'LOOKUP_TABLE default\n');
     135                                        s='%e\n';
     136                                        fprintf(fid,s,sol_struct{j}(timestep).(fieldnames{k}));
     137                                end             
     138                        end
    124139          end
    125140  end
Note: See TracChangeset for help on using the changeset viewer.