Ignore:
Timestamp:
01/30/17 16:54:51 (8 years ago)
Author:
tsantos
Message:

CHG: adjusting methods in AdaptiveMeshRefinement (not working yet).

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r21484 r21503  
    5050    if(this->fathermesh)    delete this->fathermesh;
    5151    if(this->previousmesh)  delete this->previousmesh;
    52 
     52         this->hmax=-1;
     53         this->elementswidth=-1;
    5354}
    5455/*}}}*/
    … …  
    127128/*Mesh refinement methods*/
    128129#include "TPZVTKGeoMesh.h" //itapopo
    129 void AdaptiveMeshRefinement::ExecuteRefinement(int &type_process,double *vx, double *vy, double *masklevelset, long &nvertices, long &nelements, long &nsegments, double** x, double** y, double** z, long*** elements, long*** segments){/*{{{*/
    130 
    131     //ITAPOPO _assert_(this->fathermesh);
    132     //ITAPOPO _assert_(this->previousmesh);
     130void AdaptiveMeshRefinement::ExecuteRefinement(int &type_process,double *vx, double *vy, double *masklevelset, int &nvertices, int &nelements, int &nsegments, double** x, double** y, double** z, int*** elements, int*** segments){/*{{{*/
     131
     132    //itapopo _assert_(this->fathermesh);
     133    //itapopo _assert_(this->previousmesh);
    133134   
    134135    /*Calculate the position of the grounding line using previous mesh*/
    135136    std::vector<TPZVec<REAL> > GLvec;
    136137    this->CalcGroundingLinePosition(masklevelset, GLvec);
    137 
    138    
    139     //std::ofstream file1("/Users/santos/Desktop/mesh0.vtk");
    140     //TPZVTKGeoMesh::PrintGMeshVTK(this->previousmesh,file1 );
     138   
     139    std::ofstream file1("/ronne_1/home/santos/mesh0.vtk");
     140    TPZVTKGeoMesh::PrintGMeshVTK(this->previousmesh,file1 );
    141141   
    142142    /*run refinement or unrefinement process*/
    … …  
    149149   
    150150    this->RefinementProcess(newmesh,GLvec);
    151    
    152     //std::ofstream file2("/Users/santos/Desktop/mesh1.vtk");
    153     //TPZVTKGeoMesh::PrintGMeshVTK(this->previousmesh,file2 );
     151       
     152    std::ofstream file2("/ronne_1/home/santos/mesh1.vtk");
     153    TPZVTKGeoMesh::PrintGMeshVTK(this->previousmesh,file2 );
    154154   
    155155    /*Set new mesh pointer. Previous mesh just have uniform elements*/
    … …  
    162162    TPZGeoMesh *nohangingnodesmesh = new TPZGeoMesh(*newmesh);
    163163   
    164    
    165     //std::ofstream file3("/Users/santos/Desktop/mesh2.vtk");
    166     //TPZVTKGeoMesh::PrintGMeshVTK(newmesh,file3);
     164    std::ofstream file3("/ronne_1/home/santos/mesh2.vtk");
     165    TPZVTKGeoMesh::PrintGMeshVTK(newmesh,file3);
    167166   
    168167    this->RefineMeshToAvoidHangingNodes(nohangingnodesmesh);
    169168   
    170    
    171     //std::ofstream file4("/Users/santos/Desktop/mesh3.vtk");
    172     //TPZVTKGeoMesh::PrintGMeshVTK(nohangingnodesmesh,file4);
    173    
     169    std::ofstream file4("/ronne_1/home/santos/mesh3.vtk");
     170    TPZVTKGeoMesh::PrintGMeshVTK(nohangingnodesmesh,file4);
    174171   
    175172    /*Get new geometric mesh in ISSM data structure*/
    176173    this->GetMesh(nohangingnodesmesh,nvertices,nelements,nsegments,x,y,z,elements,segments);
    177 
     174         
    178175    /*Verify the new geometry*/
    179     this->CheckMesh(nvertices,nelements,nsegments,this->elementswidth,(*x),(*y),(*z),(*elements),(*segments));
    180    
     176    this->CheckMesh(nvertices,nelements,nsegments,this->elementswidth,x,y,z,elements,segments);
     177     
    181178    delete nohangingnodesmesh;
    182 
    183179}
    184180/*}}}*/
    … …  
    189185       
    190186        /*Set elements to be refined using some criteria*/
    191         std::vector<long> ElemVec; //elements without children
     187        std::vector<int> ElemVec; //elements without children
    192188        this->SetElementsToRefine(gmesh,GLvec,hlevel,ElemVec);
    193189       
    … …  
    198194}
    199195/*}}}*/
    200 void AdaptiveMeshRefinement::RefineMesh(TPZGeoMesh *gmesh, std::vector<long> &ElemVec){/*{{{*/
     196void AdaptiveMeshRefinement::RefineMesh(TPZGeoMesh *gmesh, std::vector<int> &ElemVec){/*{{{*/
    201197
    202198        /*Refine elements in ElemVec: uniform pattern refinement*/
    203         for(long i = 0; i < ElemVec.size(); i++){
     199        for(int i = 0; i < ElemVec.size(); i++){
    204200               
    205201        /*Get geometric element and verify if it has already been refined*/
    206         long index = ElemVec[i];
     202        int index = ElemVec[i];
    207203        TPZGeoEl * geoel = gmesh->Element(index);
    208204        if(geoel->HasSubElement()) DebugStop();                              //itapopo _assert_(!geoel->HasSubElement());
    … …  
    211207        /*Divide geoel*/
    212208        TPZVec<TPZGeoEl *> Sons;
    213                 geoel->Divide(Sons);
     209                  geoel->Divide(Sons);
    214210       
    215211        /*If a 1D segment is neighbor, it must be divided too*/
    … …  
    236232   
    237233    /*Refine elements to avoid hanging nodes: non-uniform refinement*/
    238         const long NElem = gmesh->NElements();
    239     for(long i = 0; i < NElem; i++){
     234         const int NElem = gmesh->NElements();
     235    for(int i = 0; i < NElem; i++){
    240236       
    241237        /*Get geometric element and verify if it has already been refined. Geoel may not have been previously refined*/
    … …  
    253249        }
    254250       
    255         }
     251    }
    256252   
    257253    gmesh->BuildConnectivity();
    … …  
    259255}
    260256/*}}}*/
    261 void AdaptiveMeshRefinement::GetMesh(TPZGeoMesh *gmesh, long &nvertices, long &nelements, long &nsegments, double** meshX, double** meshY, double** meshZ, long*** elements, long*** segments){/*{{{*/
     257void AdaptiveMeshRefinement::GetMesh(TPZGeoMesh *gmesh, int &nvertices, int &nelements, int &nsegments, double** meshX, double** meshY, double** meshZ, int*** elements, int*** segments){/*{{{*/
    262258
    263259        /* vertices */
    264     long ntotalvertices = gmesh->NNodes();//total
     260    int ntotalvertices = gmesh->NNodes();//total
    265261   
    266262    /* mesh coords */
    … …  
    269265    double *newmeshZ = new double[ntotalvertices];
    270266       
    271         /* getting mesh coords */
    272     for(long i = 0; i < ntotalvertices; i++ ){
     267   /* getting mesh coords */
     268    for(int i = 0; i < ntotalvertices; i++ ){
    273269        TPZVec<REAL> coords(3,0.);
    274270        gmesh->NodeVec()[i].GetCoordinates(coords);
    … …  
    280276        /* elements */
    281277    std::vector<TPZGeoEl*> GeoVec; GeoVec.clear();
    282     for(long i = 0; i < gmesh->NElements(); i++){
    283    
     278    for(int i = 0; i < gmesh->NElements(); i++){
    284279        if( gmesh->ElementVec()[i]->HasSubElement() ) continue;
    285280        if( gmesh->ElementVec()[i]->MaterialId() != this->GetElemMaterialID() ) continue;
    286281        GeoVec.push_back( gmesh->ElementVec()[i]);
    287                
    288     }
    289    
    290     long ntotalelements = GeoVec.size();
    291     long ** newelements = new long*[ntotalelements];
     282    }
     283   
     284    int ntotalelements = (int)GeoVec.size();
     285    int ** newelements = new int*[ntotalelements];
    292286
    293287    if ( !(this->elementswidth == 3) && !(this->elementswidth == 4) && !(this->elementswidth == 6) ) DebugStop();
    294288
    295     for(long i = 0; i < GeoVec.size(); i++){
    296 
    297         newelements[i] = new long[this->elementswidth];
    298         for(int j = 0; j < this->elementswidth; j++) newelements[i][j] = GeoVec[i]->NodeIndex(j);
     289    for(int i = 0; i < GeoVec.size(); i++){
     290        newelements[i] = new int[this->elementswidth];
     291        for(int j = 0; j < this->elementswidth; j++) newelements[i][j] = (int)GeoVec[i]->NodeIndex(j);
    299292    }
    300293   
    301294    /* segments */
    302295    std::vector<TPZGeoEl*> SegVec; SegVec.clear();
    303     for(long i = 0; i < gmesh->NElements(); i++){
    304        
     296    for(int i = 0; i < gmesh->NElements(); i++){
    305297        if( gmesh->ElementVec()[i]->HasSubElement() ) continue;
    306298        if( gmesh->ElementVec()[i]->MaterialId() != this->GetBoundaryMaterialID() ) continue;
    307299        SegVec.push_back( gmesh->ElementVec()[i]);
    308        
    309     }
    310    
    311     long ntotalsegments = SegVec.size();
    312     long ** newsegments = new long*[ntotalsegments];
    313    
    314     for(long i = 0; i < SegVec.size(); i++){
    315        
    316         newsegments[i] = new long[3];
     300    }
     301   
     302    int ntotalsegments = (int)SegVec.size();
     303    int ** newsegments = new int*[ntotalsegments];
     304   
     305    for(int i = 0; i < SegVec.size(); i++){
     306       
     307        newsegments[i] = new int[3];
    317308        for(int j = 0; j < 2; j++) newsegments[i][j] = SegVec[i]->NodeIndex(j);
    318309       
    319         long neighborindex = SegVec[i]->Neighbour(2).Element()->Index();
    320         long neighbourid = -1;
    321        
    322         for(long j = 0; j < GeoVec.size(); j++){
     310        int neighborindex = SegVec[i]->Neighbour(2).Element()->Index();
     311        int neighbourid = -1;
     312       
     313        for(int j = 0; j < GeoVec.size(); j++){
    323314            if( GeoVec[j]->Index() == neighborindex || GeoVec[j]->FatherIndex() == neighborindex){
    324315                neighbourid = j;
    … …  
    347338    /* Find grounding line using elments center point */
    348339    GLvec.clear();
    349     for(long i=0;i<this->previousmesh->NElements();i++){
     340    for(int i=0;i<this->previousmesh->NElements();i++){
    350341       
    351342        if(this->previousmesh->Element(i)->MaterialId()!=this->GetElemMaterialID()) continue;
    … …  
    353344       
    354345        //itapopo apenas malha 2D triangular!
    355         long vertex0 = this->previousmesh->Element(i)->NodeIndex(0);
    356         long vertex1 = this->previousmesh->Element(i)->NodeIndex(1);
    357         long vertex2 = this->previousmesh->Element(i)->NodeIndex(2);
     346        int vertex0 = this->previousmesh->Element(i)->NodeIndex(0);
     347        int vertex1 = this->previousmesh->Element(i)->NodeIndex(1);
     348        int vertex2 = this->previousmesh->Element(i)->NodeIndex(2);
    358349       
    359350        double mls0 = masklevelset[vertex0];
    … …  
    384375}
    385376/*}}}*/
    386 void AdaptiveMeshRefinement::SetElementsToRefine(TPZGeoMesh *gmesh,std::vector<TPZVec<REAL> > &GLvec,int &hlevel,std::vector<long> &ElemVec){/*{{{*/
     377void AdaptiveMeshRefinement::SetElementsToRefine(TPZGeoMesh *gmesh,std::vector<TPZVec<REAL> > &GLvec,int &hlevel,std::vector<int> &ElemVec){/*{{{*/
    387378
    388379    if(!gmesh) DebugStop(); //itapopo verificar se usará _assert_
    … …  
    398389}
    399390/*}}}*/
    400 void AdaptiveMeshRefinement::TagAllElements(TPZGeoMesh *gmesh,std::vector<long> &ElemVec){/*{{{*/
     391void AdaptiveMeshRefinement::TagAllElements(TPZGeoMesh *gmesh,std::vector<int> &ElemVec){/*{{{*/
    401392   
    402393    /* Uniform refinement. This refines the entire mesh */
    403     long nelements = gmesh->NElements();
    404     for(long i=0;i<nelements;i++){
     394    int nelements = gmesh->NElements();
     395    for(int i=0;i<nelements;i++){
    405396        if(gmesh->Element(i)->MaterialId()!=this->GetElemMaterialID()) continue;
    406397        if(gmesh->Element(i)->HasSubElement()) continue;
    … …  
    409400}
    410401/*}}}*/
    411 void AdaptiveMeshRefinement::TagElementsNearGroundingLine(TPZGeoMesh *gmesh,std::vector<TPZVec<REAL> > &GLvec,int &hlevel,std::vector<long> &ElemVec){/*{{{*/
     402void AdaptiveMeshRefinement::TagElementsNearGroundingLine(TPZGeoMesh *gmesh,std::vector<TPZVec<REAL> > &GLvec,int &hlevel,std::vector<int> &ElemVec){/*{{{*/
    412403   
    413404    /* Tag elements near grounding line */
    … …  
    416407    double MaxDistance = MaxRegion / std::exp(alpha*(hlevel-1));
    417408   
    418     for(long i=0;i<gmesh->NElements();i++){
     409    for(int i=0;i<gmesh->NElements();i++){
    419410       
    420411        if(gmesh->Element(i)->MaterialId()!=this->GetElemMaterialID()) continue;
    … …  
    430421        REAL distance = MaxDistance;
    431422       
    432         for (long j = 0; j < GLvec.size(); j++) {
     423        for (int j = 0; j < GLvec.size(); j++) {
    433424           
    434425            REAL value = ( GLvec[j][0] - centerPoint[0] ) * ( GLvec[j][0] - centerPoint[0] ); // (x2-x1)^2
    … …  
    446437}
    447438/*}}}*/
    448 void AdaptiveMeshRefinement::CreateInitialMesh(long &nvertices, long &nelements, long &nsegments, int &width, double* x, double* y, double* z, long** elements, long** segments){/*{{{*/
     439void AdaptiveMeshRefinement::CreateInitialMesh(int &nvertices, int &nelements, int &nsegments, int &width, double* x, double* y, double* z, int** elements, int** segments){/*{{{*/
    449440
    450441        // itapopo _assert_(nvertices>0);
    451442    // itapoo _assert_(nelements>0);
    452     this->SetElementWidth(width);
     443   this->SetElementWidth(width);
    453444
    454445    /*Verify and creating initial mesh*/
    455446    //itapopo if(this->fathermesh) _error_("Initial mesh already exists!");
    456447   
    457     this->fathermesh = new TPZGeoMesh();
     448   this->fathermesh = new TPZGeoMesh();
    458449        this->fathermesh->NodeVec().Resize( nvertices );
    459450
    460451    /*Set the vertices (geometric nodes in NeoPZ context)*/
    461         for(long i=0;i<nvertices;i++){
     452        for(int i=0;i<nvertices;i++){
    462453       
    463454        /*x,y,z coords*/
    … …  
    469460        /*Insert in the mesh*/
    470461        this->fathermesh->NodeVec()[i].SetCoord(coord);
    471                 this->fathermesh->NodeVec()[i].SetNodeId(i);
     462                  this->fathermesh->NodeVec()[i].SetNodeId(i);
    472463        }
    473464       
    … …  
    477468    TPZManVector<long> elem(this->elementswidth,0);
    478469   
    479         for(long iel=0;iel<nelements;iel++){
     470        for(int iel=0;iel<nelements;iel++){
    480471
    481472                for(int jel=0;jel<this->elementswidth;jel++) elem[jel]=elements[iel][jel];
    … …  
    499490    TPZManVector<long> boundary(2,0.);
    500491   
    501     for(long iel=nelements;iel<nelements+nsegments;iel++){
     492    for(int iel=nelements;iel<nelements+nsegments;iel++){
    502493       
    503494        boundary[0] = segments[iel-nelements][0];
    … …  
    526517}
    527518/*}}}*/
    528 void AdaptiveMeshRefinement::CheckMesh(long &nvertices, long &nelements, long &nsegments,int &width, double* x, double* y, double* z, long** elements, long** segments){/*{{{*/
     519void AdaptiveMeshRefinement::CheckMesh(int &nvertices, int &nelements, int &nsegments,int &width, double** x, double** y, double** z, int*** elements, int*** segments){/*{{{*/
    529520
    530521    /*Basic verification*/
    … …  
    536527   
    537528    /*Verify if there are orphan nodes*/
    538     std::set<long> elemvertices;
    539     elemvertices.clear();
    540    
    541     for(long i = 0; i < nelements; i++){
    542         for(long j = 0; j < width; j++) {
    543             elemvertices.insert(elements[i][j]);
    544         }
    545     }
     529    std::set<int> elemvertices;
     530    elemvertices.clear();
     531    for(int i = 0; i < nelements; i++){
     532        for(int j = 0; j < width; j++) {
     533            elemvertices.insert((*elements)[i][j]);
     534                  }
     535         }
    546536   
    547537    if( elemvertices.size() != nvertices ) DebugStop();//itapopo verificar se irá usar o _assert_
    548    
     538       
    549539    //Verify if there are inf or NaN in coords
    550     for(long i = 0; i < nvertices; i++) if(isnan(x[i]) || isinf(x[i])) DebugStop();
    551     for(long i = 0; i < nvertices; i++) if(isnan(y[i]) || isinf(y[i])) DebugStop();
    552     for(long i = 0; i < nvertices; i++) if(isnan(z[i]) || isinf(z[i])) DebugStop();
    553     for(long i = 0; i < nvertices; i++){
    554         for(long j = 0; j < width; j++){
    555             if( isnan(elements[i][j]) || isinf(elements[i][j]) ) DebugStop();
    556         }
    557     }
    558    
    559 }
    560 /*}}}*/
     540    for(int i = 0; i < nvertices; i++) if(isnan((*x)[i]) || isinf((*x)[i])) DebugStop();
     541    for(int i = 0; i < nvertices; i++) if(isnan((*y)[i]) || isinf((*y)[i])) DebugStop();
     542    for(int i = 0; i < nvertices; i++) if(isnan((*z)[i]) || isinf((*z)[i])) DebugStop();
     543   
     544         for(int i = 0; i < nelements; i++){
     545        for(int j = 0; j < width; j++){
     546            if( isnan((*elements)[i][j]) || isinf((*elements)[i][j]) ) DebugStop();
     547        }
     548    }
     549   
     550}
     551/*}}}*/
Note: See TracChangeset for help on using the changeset viewer.