Changeset 2771


Ignore:
Timestamp:
01/06/10 09:23:45 (15 years ago)
Author:
Mathieu Morlighem
Message:

some improvements in Bamg interface (to be continued)

Location:
issm/trunk/src
Files:
4 added
6 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/Bamgx/Bamgx.cpp

    r2740 r2771  
    1111#include "../include/macros.h"
    1212#include "../toolkits/toolkits.h"
    13 #include "../EnumDefinitions/EnumDefinitions.h"
    14 
    1513
    1614/*Bamg: */
    … …  
    3028using namespace bamg;
    3129using namespace std;
    32 #define   NBVMAX 50000   
    33 
    34 
    35 int Bamgx(double** pelementsout,double** pxout,double** pyout, int* pnelout,int* pnodsout,
    36                 double* elements, double* x, double* y, int nel, int nods, double* metric, double gradation, int splitpbedge, int nbv){
    37        
    38 
     30
     31int Bamgx(BamgMesh* bamgmesh,BamgGeom* bamggeom,BamgArgs* bamgargs){
    3932       
    4033        int noerr=1;
    … …  
    4740
    4841        /*Bamg variables: */
    49         double time0, time1,time2,time3,time4;
    50         //   2 way for uses ---
    51         //   1 first mesh
    52         //   or adaptation
    53        
    54         //  double *sol=0;
    55         //  int adaptation = 0 ;
    56         Int4 i;
     42        int i,j,num;
    5743        hinterpole=1; // by def interpolation a h
    58         //  int Err=0;
    5944        int fileout=0;
    60         int nbvx = NBVMAX;
     45        int nbvx = 50000;
    6146        int iso =0,AbsError=0,nbjacoby=1,allquad=0;
    6247        int NoMeshReconstruction=0;
    … …  
    8166        double power=1;
    8267        Triangles *Thr=0,*Thb=0;
    83 
    84         char * fgeom=0,*fmeshback=0,*fmeshout=0,*fmeshr=0,*fmetrix=0,*famfmt=0,*fmsh=0,
     68        char *fmeshback=0,*fmeshout=0,*fmeshr=0,*fmetrix=0,*famfmt=0,*fmsh=0,
    8569                 *fnopo=0,*fftq=0,*fam=0,*famdba=0,*rbb=0,*rBB=0,*wbb=0,*wBB=0,
    8670                 *fMbb=0,*foM=0,*fMBB=0;
    87        
    8871        long int verbosity =10;
    89 
    9072        char *datargv[128] ;
    9173        int datargc=1;
    9274        datargv[0]= datargv[1]=0;// for create a error if no parameter
    9375
    94         //Recover command line options: */
    95         nbvx=nbv;
     76        /*testing*/
     77        int splitpbedge=1;
     78        int nel=0;
     79        int nods=0;
     80        double* x=NULL;
     81        double* y=NULL;
     82        double* elements=NULL;
     83
     84        printf("NumVertices = %i\nNumEdges = %i\n",bamggeom->NumVertices,bamggeom->NumEdges);
     85
     86        /*Recover options from inputs: */
    9687        if(splitpbedge)SplitEdgeWith2Boundary=1;
    97 
    98         /*Rest of the code taken directly from bamg.cpp, with our own mods: */
    9988
    10089        // some verification
    … …  
    121110                rBBeqMBB = !strcmp(rBB,fMBB);
    122111
    123         time0 = CPUtime();
    124 
     112        if(1){
     113                if (verbosity)
     114                 printf("Construction of a mesh from a given geometry\n");
     115
     116                Geometry Gh(bamgargs->geomfile);
     117                hmin = Max(hmin,Gh.MinimalHmin());
     118                hmax = Min(hmax,Gh.MaximalHmax());
     119               
     120                //build metric if not given in input
     121                for(i=0;i<Gh.nbv;i++)
     122                  {
     123                        MetricAnIso M=Gh[i];
     124                        MatVVP2x2 Vp(M/coef);
     125                        Vp.Maxh(hmin);
     126                        Vp.Minh(hmax);
     127                        Gh.vertices[i].m = Vp;
     128                  }
     129
     130                //generate mesh
     131                Triangles Th(nbvx,Gh);
     132
     133                /*Build output{{{2*/
     134                nelout=Th.nbt-Th.NbOutT; //number of triangles - number of external triangles
     135                nodsout=Th.nbv;
     136
     137                xout=(double*)xmalloc(nodsout*sizeof(double));
     138                yout=(double*)xmalloc(nodsout*sizeof(double));
     139                for (i=0;i<nodsout;i++){
     140                        xout[i]=Th.vertices[i].r.x;
     141                        yout[i]=Th.vertices[i].r.y;
     142                }
     143
     144                elementsout=(double*)xmalloc(3*nelout*sizeof(double));
     145                num=0;
     146                for(i=0;i<Th.nbt;i++){
     147                        Triangle &t=Th.triangles[i];
     148                        if (t.link){
     149                                //write the element only if it is part of the mesh (no boundary element)
     150                                elementsout[3*num+0]=Th.Number(t[0])+1;
     151                                elementsout[3*num+1]=Th.Number(t[1])+1;
     152                                elementsout[3*num+2]=Th.Number(t[2])+1;
     153                                num=num+1;
     154                        }
     155                }
     156
     157                /*Assign output pointers*/
     158                bamgmesh->numberofelements=nelout;
     159                bamgmesh->numberofnodes=nodsout;
     160                bamgmesh->index=elementsout;
     161                bamgmesh->x=xout;
     162                bamgmesh->y=yout;
     163
     164                return noerr;
     165                /*}}}2*/
     166        }
     167
     168        if (0){
     169        /*Anisotropic mesh adaptation {{{1*/
    125170        /*Read background mesh from simple delaunay triangulation: */
    126171        Triangles BTh(elements,x,y,nel,nods,nbvx,cutoffradian); // read the background mesh
    … …  
    132177
    133178        BTh.MakeQuadTree();
    134         time1 = CPUtime();
    135179        if (fmetrix)
    136180                BTh.ReadMetric(fmetrix,hmin,hmax,coef);
    … …  
    182226                anisomax=1;
    183227        BTh.BoundAnisotropy(anisomax,hminaniso);
    184         time2 = CPUtime();
    185228        if(foM)
    186229        {
    … …  
    215258                        Th.SplitInternalEdgeWithBorderVertices();
    216259                Th.ReNumberingTheTriangleBySubDomain();
    217                 time3 = CPUtime();
    218 
    219                 if(verbosity>0){
    220                         cout << " Cpu time for meshing          : " << time3-time2 << "s  Nb Vertices/s = " <<  (Th.nbv) /(time3-time2) << endl;
    221                 }
     260
    222261                if(verbosity>3) Th.ShowHistogram();
    223262                if(NbSmooth>0) Th.SmoothingVertex(NbSmooth,omega);
    … …  
    243282                                        cout << " -- interpolation P1  files : " << rBB<< " -> " << wBB <<endl;
    244283                        }
    245                         double time00  = CPUtime();       
    246284                        const int dim=2;
    247285                        // optimisation read only si rbb != fMbb
    … …  
    336374                                if (solBB)
    337375                                        delete [] solBB;
    338                                 double  time04 = CPUtime();
    339                                 if(verbosity>0)
    340                                         cout << " Cpu time for interpolation " << time04-time00 << " s" <<endl;
    341 
    342                 }
    343                 time4 = CPUtime();
     376
     377                }
    344378
    345379                if(verbosity>0)
    346380                {
    347                         cout << " Cpu time for meshing with io  : " << time4-time0
    348                                 << "s  Nb Vertices/s = " << Th.nbv/(time4-time0)
    349                                 << endl
    350                                 << " Nb vertices = " << Th.nbv;
    351381                        if (Th.nbt-Th.NbOutT-Th.NbOfQuad*2)
    352382                                cout  << " Nb Triangles = " << (Th.nbt-Th.NbOutT-Th.NbOfQuad*2);
    … …  
    356386                }
    357387
    358 
    359388                // cout << "delete &Th " ;
    360389                delete &Th;
    361390                //cout << " end " << endl;
    362391        }
    363 
    364         throw ErrorException(__FUNCT__," almost there\n");
    365 
    366 
    367392        /*Assign output pointers:*/
    368         *pelementsout=elementsout;
    369         *pxout=xout;
    370         *pyout=yout;
    371         *pnelout=nelout;
    372         *pnodsout=nodsout;
     393        bamgmesh->numberofelements=nelout;
     394        bamgmesh->numberofnodes=nodsout;
     395        bamgmesh->index=elementsout;
     396        bamgmesh->x=xout;
     397        bamgmesh->y=yout;
    373398
    374399        return noerr;
     400        /*}}}*/
     401        }
    375402}

  • issm/trunk/src/c/Bamgx/Bamgx.h

    r2740 r2771  
    99
    1010/* local prototypes: */
    11 int     Bamgx(double** pelementsout,double** pxout,double** pyout, int* pnelout,int* pnodsout,
    12                 double* elements, double* x, double* y, int nel, int nods, double* metric, double gradation, int splitpbedge, int nbv);
     11int     Bamgx(BamgMesh* bamgmesh,BamgGeom* bamggeom,BamgArgs* bamgargs);
    1312
    1413#endif  /* _BAMGX_H */

  • issm/trunk/src/c/Bamgx/MeshRead.cpp

    r2740 r2771  
    1 // -*- Mode : c++ -*-
    2 //
    3 // SUMMARY  :     
    4 // USAGE    :       
    5 // ORG      :
    6 // AUTHOR   : Frederic Hecht
    7 // E-MAIL   : hecht@ann.jussieu.fr
    8 //
    9 
    10 /*
    11  
    12  This file is part of Freefem++
     1/*This file is largely inspired by Freefem++
    132 
    143 Freefem++ is free software; you can redistribute it and/or modify

  • issm/trunk/src/c/objects/objects.h

    r1843 r2771  
    3939#include "./OptArgs.h"
    4040#include "./OptPars.h"
     41#include "./BamgGeom.h"
     42#include "./BamgMesh.h"
     43#include "./BamgArgs.h"
    4144
    4245

  • issm/trunk/src/mex/Bamg/Bamg.cpp

    r2739 r2771  
    22 *\brief: bamg module.
    33 */
    4 
    54#include "./Bamg.h"
    6 
    75
    86void mexFunction( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]){
    … …  
    108        /*diverse: */
    119        int   noerr=1;
     10        int   i;
     11        BamgArgs bamgargs;
     12        BamgMesh bamgmesh;
     13        BamgGeom bamggeom;
    1214
    1315        /*inputs: */
    14         double* elements=NULL;
    15         int     dummy;
    16         double* x=NULL;
    17         double* y=NULL;
    18         double* metric=NULL;
    19         int     nods,nel;
    20         double  gradation;
    21         int     splitpbedge;
    22         int     nbv;
     16        char*   geomfile=NULL;
     17        int     NumVertices;
     18        double* Vertices=NULL;
     19        int     NumEdges;
     20        double* Edges=NULL;
     21        double* hVertices=NULL;
    2322
    24         /*output: */
    25         double* elementsout=NULL;
    26         double* xout=NULL;
    27         double* yout=NULL;
    28         int     nodsout,nelout;
    2923
    3024        /*Boot module: */
    … …  
    3428        CheckNumMatlabArguments(nlhs,NLHS,nrhs,NRHS,__FUNCT__,&BamgUsage);
    3529
    36         /*Input datasets: */
    37         FetchData(&elements,&nel,&dummy,ELEMENTS);
    38         FetchData(&x,&nods,X);
    39         FetchData(&y,&nods,Y);
    40         FetchData(&metric,&nods,METRIC);
    41         FetchData(&gradation,GRADATION);
    42         FetchData(&splitpbedge,SPLITPBEDGE);
    43         FetchData(&nbv,NBV);
    44        
     30        /*process inputs*/
     31        FetchData(&geomfile,mxGetField(BAMGOPTIONS,0,"fgeom"));
     32
     33        FetchData(&NumVertices,mxGetField(BAMGGEOMETRY,0,"NumVertices"));
     34        FetchData(&Vertices,NULL,NULL,mxGetField(BAMGGEOMETRY,0,"Vertices"));
     35        FetchData(&NumEdges,mxGetField(BAMGGEOMETRY,0,"NumEdges"));
     36        FetchData(&Edges,NULL,NULL,mxGetField(BAMGGEOMETRY,0,"Edges"));
     37        FetchData(&hVertices,NULL,NULL,mxGetField(BAMGGEOMETRY,0,"hVertices"));
     38
     39        /*create bamg geometry input*/
     40        bamgargs.geomfile=geomfile;
     41
     42        bamggeom.NumVertices=NumVertices;
     43        bamggeom.Vertices=Vertices;
     44        bamggeom.NumEdges=NumEdges;
     45        bamggeom.Edges=Edges;
     46        bamggeom.hVertices=hVertices;
     47
    4548        /*!Generate internal degree of freedom numbers: */
    46         Bamgx(&elementsout,&xout,&yout,&nelout,&nodsout,elements,x,y,nel,nods,metric,gradation,splitpbedge,nbv);
     49        Bamgx(&bamgmesh,&bamggeom,&bamgargs);
    4750
    4851        /*write output datasets: */
    49         WriteData(ELEMENTSOUT,elementsout,nelout,3);
    50         WriteData(XOUT,xout,nodsout);
    51         WriteData(YOUT,yout,nodsout);
     52        WriteData(ELEMENTSOUT,bamgmesh.index,bamgmesh.numberofelements,3);
     53        WriteData(XOUT,bamgmesh.x,bamgmesh.numberofnodes);
     54        WriteData(YOUT,bamgmesh.y,bamgmesh.numberofnodes);
    5255
    5356        /*Free ressources: */
    54         xfree((void**)&elements);
    55         xfree((void**)&x);
    56         xfree((void**)&y);
    57         xfree((void**)&metric);
    58         xfree((void**)&elementsout);
    59         xfree((void**)&xout);
    60         xfree((void**)&yout);
     57        xfree((void**)&geomfile);
     58        xfree((void**)&Vertices);
     59        xfree((void**)&Edges);
     60        xfree((void**)&hVertices);
    6161
    6262        /*end module: */

  • issm/trunk/src/mex/Bamg/Bamg.h

    r2739 r2771  
    1818   
    1919/* serial input macros: */
    20 #define ELEMENTS (mxArray*)prhs[0]
    21 #define X (mxArray*)prhs[1]
    22 #define Y (mxArray*)prhs[2]
    23 #define METRIC (mxArray*)prhs[3]
    24 #define GRADATION (mxArray*)prhs[4]
    25 #define SPLITPBEDGE (mxArray*)prhs[5]
    26 #define NBV (mxArray*)prhs[6]
     20#define BAMGGEOMETRY (mxArray*)prhs[0]
     21#define BAMGOPTIONS  (mxArray*)prhs[1]
    2722
    2823/* serial output macros: */
    … …  
    3530#define NLHS  3
    3631#undef NRHS
    37 #define NRHS  7
     32#define NRHS  2
    3833
    3934#endif  /* _BAMG_H */
Note: See TracChangeset for help on using the changeset viewer.