Changeset 2780


Ignore:
Timestamp:
01/07/10 09:31:16 (15 years ago)
Author:
Mathieu Morlighem
Message:

some improvements in Bamg

Location:
issm/trunk/src
Files:
5 edited

Legend:

Unmodified
Added
Removed

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

    r2772 r2780  
    2020#include <new>
    2121#include <cassert>
     22#include <iomanip>
     23#include <fstream>
    2224#include "Meshio.h"
    23 #include <iomanip>
    2425#include "Mesh2.h"
    2526#include "QuadTree.h"
    26 #include <fstream>
    2727
    2828using namespace bamg;
    … …  
    4343        hinterpole=1; // by def interpolation a h
    4444        int fileout=0;
    45         int nbvx = 50000;
     45        int nbvx = 1000000;
    4646        int iso =0,AbsError=0,nbjacoby=1,allquad=0;
    4747        int NoMeshReconstruction=0;
    … …  
    5050        Real8 cutoffradian=-1;
    5151        double anisomax = 1e6;
    52         double err=0.01,errg=0.1,coef=1,hmin=1.e-100,hmax=1.e17,raison=0,CutOff=1e-5;
     52        double err=0.01,errg=0.1,coef=1,hmin=1.e-100,hmax=1.e17,gradation=1.2,cutoff=1e-5;
    5353        int KeepBackVertices=1;
    5454        double hminaniso=1e-100;
    … …  
    7676        /*testing*/
    7777        int splitpbedge=1;
    78         int nel=0;
    79         int nods=0;
    80         double* x=NULL;
    81         double* y=NULL;
    82         double* elements=NULL;
     78
     79        /*Bamg options*/
     80        iso=bamgargs->iso;
     81        hmin=bamgargs->hmin;
     82        hmax=bamgargs->hmax;
     83        gradation=bamgargs->gradation;
     84        cutoff=bamgargs->cutoff;
    8385
    8486        /*Recover options from inputs: */
    … …  
    108110                rBBeqMBB = !strcmp(rBB,fMBB);
    109111
    110         if(1){
     112        if(bamgmesh->numberofelements==0){
     113                /*Mesh generation {{{1*/
    111114                if (verbosity>0) printf("Construction of a mesh from a given geometry\n");
    112115                Geometry Gh(bamggeom);
    113116                hmin = Max(hmin,Gh.MinimalHmin());
    114117                hmax = Min(hmax,Gh.MaximalHmax());
    115                
     118
    116119                //build metric if not given in input
    117120                for(i=0;i<Gh.nbv;i++)
    … …  
    127130                Triangles Th(nbvx,Gh);
    128131
    129                 /*Build output {{{1*/
     132                //Build output
    130133                nelout=Th.nbt-Th.NbOutT; //number of triangles - number of external triangles
    131134                nodsout=Th.nbv;
    … …  
    150153                        }
    151154                }
    152 
    153                 /*Assign output pointers*/
    154                 bamgmesh->numberofelements=nelout;
    155                 bamgmesh->numberofnodes=nodsout;
    156                 bamgmesh->index=elementsout;
    157                 bamgmesh->x=xout;
    158                 bamgmesh->y=yout;
    159 
    160                 return noerr;
    161                 /*}}}1*/
    162 
     155                /*}}}*/
    163156        }
    164157
    165         if (0){
    166         /*Anisotropic mesh adaptation {{{1*/
    167         /*Read background mesh from simple delaunay triangulation: */
    168         Triangles BTh(elements,x,y,nel,nods,nbvx,cutoffradian); // read the background mesh
    169 
    170         hmin = Max(hmin,BTh.MinimalHmin());
    171         hmax = Min(hmax,BTh.MaximalHmax());
    172 
    173         throw ErrorException(__FUNCT__," done");
    174 
    175         BTh.MakeQuadTree();
    176         if (fmetrix)
    177                 BTh.ReadMetric(fmetrix,hmin,hmax,coef);
    178         else
    179         { // init with hmax
    180                 Metric Mhmax(hmax);
    181                 for (Int4 iv=0;iv<BTh.nbv;iv++)
    182                         BTh[iv].m = Mhmax;
     158        else{
     159                /*Anisotropic mesh adaptation {{{1*/
     160                if (verbosity>0) printf("Anisotropic mesh adaptation\n");
     161
     162                /*Read background mesh from simple delaunay triangulation: */
     163                Triangles BTh(bamgmesh->index,bamgmesh->x,bamgmesh->y,bamgmesh->numberofelements,bamgmesh->numberofnodes,nbvx,cutoffradian); // read the background mesh
     164
     165                hmin = Max(hmin,BTh.MinimalHmin());
     166                hmax = Min(hmax,BTh.MaximalHmax());
     167
     168                throw ErrorException(__FUNCT__," done");
     169
     170                BTh.MakeQuadTree();
     171                if (fmetrix){
     172                        BTh.ReadMetric(fmetrix,hmin,hmax,coef);
     173                }
     174                else { // init with hmax
     175                        Metric Mhmax(hmax);
     176                        for (Int4 iv=0;iv<BTh.nbv;iv++)
     177                         BTh[iv].m = Mhmax;
     178                }
     179                if (fMbb) {
     180                        solMbb = ReadbbFile(fMbb,nbsolbb,lsolbb,2,2);
     181                        if (lsolbb != BTh.nbv)
     182                          {
     183                                cerr << "fatal error  nbsol " << nbsolbb << " " << lsolbb<< " =! " << BTh.nbv << endl;
     184                                cerr << "  size of sol incorrect " << endl;
     185                                MeshError(99);
     186                          }
     187                        assert(lsolbb==BTh.nbv);
     188                        BTh.IntersectConsMetric(solMbb,nbsolbb,0,hmin,hmax,sqrt(err)*coef,1e30,AbsError?0.0:cutoff,
     189                                                nbjacoby,Rescaling,power,ChoiseHessien);
     190                        if(!rbbeqMbb)
     191                         delete [] solMbb,solMbb =0;
     192
     193                }
     194                if (fMBB) {
     195                        solMBB = ReadBBFile(fMBB,nbsolBB,lsolBB,typesolsBB,2,2);
     196                        if (lsolBB != BTh.nbv)
     197                          {
     198                                cerr << "fatal error  nbsol " << nbsolBB << " " << lsolBB << " =! " << BTh.nbv << endl;
     199                                cerr << "  size of sol incorrect " << endl;
     200                                MeshError(99);
     201                          }
     202                        assert(lsolBB==BTh.nbv);
     203                        BTh.IntersectConsMetric(solMBB,nbsolBB,0,hmin,hmax,sqrt(err)*coef,1e30,AbsError?0.0:cutoff,
     204                                                nbjacoby,Rescaling,ChoiseHessien);
     205                        if(!rBBeqMBB)
     206                         delete [] solMBB,solMBB =0;
     207
     208                }
     209
     210                BTh.IntersectGeomMetric(errg,iso);
     211
     212                if(gradation) BTh.SmoothMetric(gradation);
     213                BTh.MaxSubDivision(maxsubdiv);
     214                BTh.BoundAnisotropy(anisomax,hminaniso);
     215
     216                if(foM) {
     217                        if(verbosity >2)
     218                         cout << " -- write Metric  file " << foM <<endl;
     219
     220                        ofstream f(foM);
     221                        if(f) BTh.WriteMetric(f,iso);
     222                }
     223
     224                if ( fileout) {
     225                        if ( NoMeshReconstruction)
     226                         if (( fmeshback == fmeshr) || (!strcmp(fmeshback,fmeshr)))
     227                          Thr=&BTh,Thb=0; // back and r mesh are the same
     228                         else
     229                          Thr= new Triangles(fmeshr,cutoffradian),Thb=&BTh; // read the new
     230
     231
     232
     233                        Triangles & Th( *(NoMeshReconstruction
     234                                                        ?  new Triangles(*Thr,&Thr->Gh,Thb,nbvx) // copy the mesh + free space to modification 
     235                                                        :  new Triangles(nbvx,BTh,KeepBackVertices)     // construct a new mesh
     236                                                        ));
     237                        if (Thr != &BTh) delete Thr;
     238
     239                        if(costheta<=1.0)
     240                         Th.MakeQuadrangles(costheta);
     241                        if (allquad)
     242                         Th.SplitElement(allquad==2);
     243                        if(SplitEdgeWith2Boundary)
     244                         Th.SplitInternalEdgeWithBorderVertices();
     245                        Th.ReNumberingTheTriangleBySubDomain();
     246
     247                        if(verbosity>3) Th.ShowHistogram();
     248                        if(NbSmooth>0) Th.SmoothingVertex(NbSmooth,omega);
     249                        if(verbosity>3 && NbSmooth>0) Th.ShowHistogram();
     250
     251                        Th.BTh.ReMakeTriangleContainingTheVertex();
     252
     253                        if (fmeshout) Th.Write(fmeshout  ,Triangles::BDMesh);
     254                        if (famfmt)   Th.Write(famfmt    ,Triangles::am_fmtMesh);
     255                        if (fam)      Th.Write(fam       ,Triangles::amMesh);
     256                        if (famdba)   Th.Write(famdba    ,Triangles::amdbaMesh);
     257                        if (fftq)     Th.Write(fftq      ,Triangles::ftqMesh);
     258                        if (fmsh)     Th.Write(fmsh      ,Triangles::mshMesh);
     259                        if (fnopo)    Th.Write(fnopo     ,Triangles::NOPOMesh);
     260
     261                        if ( ( rbb && wbb)  ||( rBB && wBB)){  // the code for interpolation
     262
     263                                if(verbosity >1)
     264                                  {
     265                                        if (rbb )
     266                                         cout << " -- interpolation P1  files : " << rbb << " -> " << wbb <<endl;
     267                                        if (rBB )
     268                                         cout << " -- interpolation P1  files : " << rBB<< " -> " << wBB <<endl;
     269                                  }
     270                                const int dim=2;
     271                                // optimisation read only si rbb != fMbb
     272
     273                                double *solbb=0;
     274                                double *solBB=0;
     275
     276                                if (rbb)
     277                                 solbb =  rbbeqMbb? solMbb : ReadbbFile(rbb,nbsolbb,lsolbb,2,2);
     278                                if (rBB)
     279                                 solBB =  rBBeqMBB? solMBB : ReadBBFile(rBB,nbsolBB,lsolBB,typesolsBB,2,2);
     280
     281                                // cout << " " << rbbeqMbb << " " <<  sol << " " << nbsol << " " << lsol << endl;
     282                                if (!solBB && !solbb )
     283                                  {
     284                                        cerr << " Fatal Error "  << "solBB =  " <<  solBB << " solbb= " << solbb << endl;
     285                                        exit(2);}
     286
     287                                        ofstream *fbb = wbb ? new ofstream(wbb) :0;
     288                                        ofstream *fBB = wBB ? new ofstream(wBB) :0;
     289                                        Int4   nbfieldBB = 0, nbfieldbb = nbsolbb;
     290                                        if (fbb)
     291                                         *fbb  << dim << " " << nbsolbb << " " << Th.nbv <<" " << 2 << endl;
     292                                        if (fBB)
     293                                          {
     294                                                int i;
     295                                                *fBB  << dim << " " << nbsolBB ;
     296                                                for ( i=0;i<nbsolBB;i++)
     297                                                 *fBB << " " << (typesolsBB ?typesolsBB[i]+1:1) ;
     298                                                *fBB << " " << Th.nbv <<" " << 2 << endl;
     299                                                assert(dim==2);
     300                                                for ( i=0;i<nbsolBB;i++)
     301                                                 nbfieldBB += typesolsBB ? typesolsBB[i]+1 : 1;
     302                                                // this code is good only if dim == 2
     303                                          }
     304                                        cout << "nb of field BB " << nbfieldBB << endl;
     305                                        //               if(fBB) fBB->precision(15);
     306                                        //if(fbb) fbb->precision(15);
     307                                        for(i=0;i<Th.nbv;i++)
     308                                          {
     309                                                Int4 i0,i1,i2;
     310                                                double a[3];
     311                                                Icoor2 dete[3];
     312                                                I2 I = Th.BTh.toI2(Th.vertices[i].r);
     313                                                Triangle & tb = *Th.BTh.FindTriangleContening(I,dete);
     314
     315                                                if (tb.det>0) { // internal point
     316                                                        a[0]= (Real8) dete[0]/ tb.det;
     317                                                        a[1]= (Real8) dete[1] / tb.det;
     318                                                        a[2] = (Real8) dete[2] / tb.det;
     319                                                        i0=Th.BTh.Number(tb[0]);
     320                                                        i1=Th.BTh.Number(tb[1]);
     321                                                        i2=Th.BTh.Number(tb[2]);}
     322                                                else {
     323                                                        double aa,bb;
     324
     325                                                        TriangleAdjacent  ta=CloseBoundaryEdge(I,&tb,aa,bb).Adj();
     326
     327                                                        int k = ta;
     328                                                        Triangle & tc = *(Triangle *)ta;
     329                                                        i0=Th.BTh.Number(tc[0]);
     330                                                        i1=Th.BTh.Number(tc[1]);
     331                                                        i2=Th.BTh.Number(tc[2]);
     332                                                        a[VerticesOfTriangularEdge[k][1]] =aa;
     333                                                        a[VerticesOfTriangularEdge[k][0]] = bb;
     334                                                        a[OppositeVertex[k]] = 1- aa -bb;}
     335
     336                                                        Int4  ibb0 = nbfieldbb*i0;
     337                                                        Int4  ibb1 = nbfieldbb*i1;
     338                                                        Int4  ibb2 = nbfieldbb*i2;
     339
     340                                                        Int4  iBB0 = nbfieldBB*i0;
     341                                                        Int4  iBB1 = nbfieldBB*i1;
     342                                                        Int4  iBB2 = nbfieldBB*i2;
     343                                                        Int4 j=0;
     344
     345                                                        for ( j=0;j<nbfieldbb;j++)
     346                                                         *fbb  << " " << ( a[0] * solbb[ibb0++] + a[1] * solbb[ibb1++] + a[2]* solbb[ibb2++]) ;
     347
     348                                                        for (j=0;j<nbfieldBB;j++)
     349                                                         *fBB  << " " << ( a[0] * solBB[iBB0++] + a[1] * solBB[iBB1++] + a[2]* solBB[iBB2++]) ;
     350
     351                                                        if (fbb) *fbb  << endl;
     352                                                        if (fBB) *fBB  << endl;
     353                                          }
     354                                        if (fbb)
     355                                         delete fbb  ; // close
     356                                        if (fBB)
     357                                         delete fBB  ; // close
     358                                        if (solbb)
     359                                         delete [] solbb;
     360                                        if (solBB)
     361                                         delete [] solBB;
     362
     363                        }
     364
     365                        if(verbosity>0)
     366                          {
     367                                if (Th.nbt-Th.NbOutT-Th.NbOfQuad*2)
     368                                 cout  << " Nb Triangles = " << (Th.nbt-Th.NbOutT-Th.NbOfQuad*2);
     369                                if (Th.NbOfQuad )
     370                                 cout  << " Nb Quadrilaterals = " << Th.NbOfQuad  ;
     371                                cout   << endl;
     372                          }
     373
     374                        //Build output
     375                        nelout=Th.nbt-Th.NbOutT; //number of triangles - number of external triangles
     376                        nodsout=Th.nbv;
     377
     378                        xout=(double*)xmalloc(nodsout*sizeof(double));
     379                        yout=(double*)xmalloc(nodsout*sizeof(double));
     380                        for (i=0;i<nodsout;i++){
     381                                xout[i]=Th.vertices[i].r.x;
     382                                yout[i]=Th.vertices[i].r.y;
     383                        }
     384
     385                        elementsout=(double*)xmalloc(3*nelout*sizeof(double));
     386                        num=0;
     387                        for(i=0;i<Th.nbt;i++){
     388                                Triangle &t=Th.triangles[i];
     389                                if (t.link){
     390                                        //write the element only if it is part of the mesh (no boundary element)
     391                                        elementsout[3*num+0]=Th.Number(t[0])+1;
     392                                        elementsout[3*num+1]=Th.Number(t[1])+1;
     393                                        elementsout[3*num+2]=Th.Number(t[2])+1;
     394                                        num=num+1;
     395                                }
     396                        }
     397
     398                        /*clean up*/
     399                        delete &Th;
     400                }
     401                /*}}}*/
    183402        }
    184         if (fMbb)
    185         {
    186                 solMbb = ReadbbFile(fMbb,nbsolbb,lsolbb,2,2);
    187                 if (lsolbb != BTh.nbv)
    188                 {
    189                         cerr << "fatal error  nbsol " << nbsolbb << " " << lsolbb<< " =! " << BTh.nbv << endl;
    190                         cerr << "  size of sol incorrect " << endl;
    191                         MeshError(99);
    192                 }
    193                 assert(lsolbb==BTh.nbv);
    194                 BTh.IntersectConsMetric(solMbb,nbsolbb,0,hmin,hmax,sqrt(err)*coef,1e30,AbsError?0.0:CutOff,
    195                                 nbjacoby,Rescaling,power,ChoiseHessien);
    196                 if(!rbbeqMbb)
    197                         delete [] solMbb,solMbb =0;
    198 
    199         }
    200         if (fMBB)
    201         {
    202                 solMBB = ReadBBFile(fMBB,nbsolBB,lsolBB,typesolsBB,2,2);
    203                 if (lsolBB != BTh.nbv)
    204                 {
    205                         cerr << "fatal error  nbsol " << nbsolBB << " " << lsolBB << " =! " << BTh.nbv << endl;
    206                         cerr << "  size of sol incorrect " << endl;
    207                         MeshError(99);
    208                 }
    209                 assert(lsolBB==BTh.nbv);
    210                 BTh.IntersectConsMetric(solMBB,nbsolBB,0,hmin,hmax,sqrt(err)*coef,1e30,AbsError?0.0:CutOff,
    211                                 nbjacoby,Rescaling,ChoiseHessien);
    212                 if(!rBBeqMBB)
    213                         delete [] solMBB,solMBB =0;
    214 
    215         }
    216 
    217         BTh.IntersectGeomMetric(errg,iso);
    218         if(raison)
    219                 BTh.SmoothMetric(raison);
    220         BTh.MaxSubDivision(maxsubdiv);
    221         //
    222         if (iso)
    223                 anisomax=1;
    224         BTh.BoundAnisotropy(anisomax,hminaniso);
    225         if(foM)
    226         {
    227                 if(verbosity >2)
    228                         cout << " -- write Metric  file " << foM <<endl;
    229 
    230                 ofstream f(foM);
    231                 if(f) BTh.WriteMetric(f,iso);
    232         }
    233 
    234         if ( fileout)
    235         {
    236                 if ( NoMeshReconstruction)
    237                         if (( fmeshback == fmeshr) || (!strcmp(fmeshback,fmeshr)))
    238                                 Thr=&BTh,Thb=0; // back and r mesh are the same
    239                         else
    240                                 Thr= new Triangles(fmeshr,cutoffradian),Thb=&BTh; // read the new
    241 
    242 
    243 
    244                 Triangles & Th( *(NoMeshReconstruction
    245                                         ?  new Triangles(*Thr,&Thr->Gh,Thb,nbvx) // copy the mesh + free space to modification 
    246                                         :  new Triangles(nbvx,BTh,KeepBackVertices)     // construct a new mesh
    247                                         ));
    248                 if (Thr != &BTh) delete Thr;
    249 
    250                 if(costheta<=1.0)
    251                         Th.MakeQuadrangles(costheta);
    252                 if (allquad)
    253                         Th.SplitElement(allquad==2);
    254                 if(SplitEdgeWith2Boundary)
    255                         Th.SplitInternalEdgeWithBorderVertices();
    256                 Th.ReNumberingTheTriangleBySubDomain();
    257 
    258                 if(verbosity>3) Th.ShowHistogram();
    259                 if(NbSmooth>0) Th.SmoothingVertex(NbSmooth,omega);
    260                 if(verbosity>3 && NbSmooth>0) Th.ShowHistogram();
    261                
    262                 Th.BTh.ReMakeTriangleContainingTheVertex();
    263 
    264                 if (fmeshout) Th.Write(fmeshout  ,Triangles::BDMesh);
    265                 if (famfmt)   Th.Write(famfmt    ,Triangles::am_fmtMesh);
    266                 if (fam)      Th.Write(fam       ,Triangles::amMesh);
    267                 if (famdba)   Th.Write(famdba    ,Triangles::amdbaMesh);
    268                 if (fftq)     Th.Write(fftq      ,Triangles::ftqMesh);
    269                 if (fmsh)     Th.Write(fmsh      ,Triangles::mshMesh);
    270                 if (fnopo)    Th.Write(fnopo     ,Triangles::NOPOMesh);
    271 
    272                 if ( ( rbb && wbb)  ||( rBB && wBB)){  // the code for interpolation
    273 
    274                         if(verbosity >1)
    275                         {
    276                                 if (rbb )
    277                                         cout << " -- interpolation P1  files : " << rbb << " -> " << wbb <<endl;
    278                                 if (rBB )
    279                                         cout << " -- interpolation P1  files : " << rBB<< " -> " << wBB <<endl;
    280                         }
    281                         const int dim=2;
    282                         // optimisation read only si rbb != fMbb
    283 
    284                         double *solbb=0;
    285                         double *solBB=0;
    286 
    287                         if (rbb)
    288                                 solbb =  rbbeqMbb? solMbb : ReadbbFile(rbb,nbsolbb,lsolbb,2,2);
    289                         if (rBB)
    290                                 solBB =  rBBeqMBB? solMBB : ReadBBFile(rBB,nbsolBB,lsolBB,typesolsBB,2,2);
    291 
    292                         // cout << " " << rbbeqMbb << " " <<  sol << " " << nbsol << " " << lsol << endl;
    293                         if (!solBB && !solbb )
    294                         {
    295                                 cerr << " Fatal Error "  << "solBB =  " <<  solBB << " solbb= " << solbb << endl;
    296                                 exit(2);}
    297 
    298                                 ofstream *fbb = wbb ? new ofstream(wbb) :0;
    299                                 ofstream *fBB = wBB ? new ofstream(wBB) :0;
    300                                 Int4   nbfieldBB = 0, nbfieldbb = nbsolbb;
    301                                 if (fbb)
    302                                         *fbb  << dim << " " << nbsolbb << " " << Th.nbv <<" " << 2 << endl;
    303                                 if (fBB)
    304                                 {
    305                                         int i;
    306                                         *fBB  << dim << " " << nbsolBB ;
    307                                         for ( i=0;i<nbsolBB;i++)
    308                                                 *fBB << " " << (typesolsBB ?typesolsBB[i]+1:1) ;
    309                                         *fBB << " " << Th.nbv <<" " << 2 << endl;
    310                                         assert(dim==2);
    311                                         for ( i=0;i<nbsolBB;i++)
    312                                                 nbfieldBB += typesolsBB ? typesolsBB[i]+1 : 1;
    313                                         // this code is good only if dim == 2
    314                                 }
    315                                 cout << "nb of field BB " << nbfieldBB << endl;
    316                                 //               if(fBB) fBB->precision(15);
    317                                 //if(fbb) fbb->precision(15);
    318                                 for(i=0;i<Th.nbv;i++)
    319                                 {
    320                                         Int4 i0,i1,i2;
    321                                         double a[3];
    322                                         Icoor2 dete[3];
    323                                         I2 I = Th.BTh.toI2(Th.vertices[i].r);
    324                                         Triangle & tb = *Th.BTh.FindTriangleContening(I,dete);
    325 
    326                                         if (tb.det>0) { // internal point
    327                                                 a[0]= (Real8) dete[0]/ tb.det;
    328                                                 a[1]= (Real8) dete[1] / tb.det;
    329                                                 a[2] = (Real8) dete[2] / tb.det;
    330                                                 i0=Th.BTh.Number(tb[0]);
    331                                                 i1=Th.BTh.Number(tb[1]);
    332                                                 i2=Th.BTh.Number(tb[2]);}
    333                                         else {
    334                                                 double aa,bb;
    335 
    336                                                 TriangleAdjacent  ta=CloseBoundaryEdge(I,&tb,aa,bb).Adj();
    337 
    338                                                 int k = ta;
    339                                                 Triangle & tc = *(Triangle *)ta;
    340                                                 i0=Th.BTh.Number(tc[0]);
    341                                                 i1=Th.BTh.Number(tc[1]);
    342                                                 i2=Th.BTh.Number(tc[2]);
    343                                                 a[VerticesOfTriangularEdge[k][1]] =aa;
    344                                                 a[VerticesOfTriangularEdge[k][0]] = bb;
    345                                                 a[OppositeVertex[k]] = 1- aa -bb;}
    346 
    347                                                 Int4  ibb0 = nbfieldbb*i0;
    348                                                 Int4  ibb1 = nbfieldbb*i1;
    349                                                 Int4  ibb2 = nbfieldbb*i2;
    350 
    351                                                 Int4  iBB0 = nbfieldBB*i0;
    352                                                 Int4  iBB1 = nbfieldBB*i1;
    353                                                 Int4  iBB2 = nbfieldBB*i2;
    354                                                 Int4 j=0;
    355 
    356                                                 for ( j=0;j<nbfieldbb;j++)
    357                                                         *fbb  << " " << ( a[0] * solbb[ibb0++] + a[1] * solbb[ibb1++] + a[2]* solbb[ibb2++]) ;
    358 
    359                                                 for (j=0;j<nbfieldBB;j++)
    360                                                         *fBB  << " " << ( a[0] * solBB[iBB0++] + a[1] * solBB[iBB1++] + a[2]* solBB[iBB2++]) ;
    361 
    362                                                 if (fbb) *fbb  << endl;
    363                                                 if (fBB) *fBB  << endl;
    364                                 }
    365                                 if (fbb)
    366                                         delete fbb  ; // close
    367                                 if (fBB)
    368                                         delete fBB  ; // close
    369                                 if (solbb)
    370                                         delete [] solbb;
    371                                 if (solBB)
    372                                         delete [] solBB;
    373 
    374                 }
    375 
    376                 if(verbosity>0)
    377                 {
    378                         if (Th.nbt-Th.NbOutT-Th.NbOfQuad*2)
    379                                 cout  << " Nb Triangles = " << (Th.nbt-Th.NbOutT-Th.NbOfQuad*2);
    380                         if (Th.NbOfQuad )
    381                                 cout  << " Nb Quadrilaterals = " << Th.NbOfQuad  ;
    382                         cout   << endl;
    383                 }
    384 
    385                 // cout << "delete &Th " ;
    386                 delete &Th;
    387                 //cout << " end " << endl;
    388         }
    389         /*Assign output pointers:*/
     403
     404        /*Assign output pointers*/
    390405        bamgmesh->numberofelements=nelout;
    391406        bamgmesh->numberofnodes=nodsout;
     407        xfree((void**)&bamgmesh->index);
    392408        bamgmesh->index=elementsout;
     409        xfree((void**)&bamgmesh->x);
    393410        bamgmesh->x=xout;
     411        xfree((void**)&bamgmesh->y);
    394412        bamgmesh->y=yout;
    395413
     414        /*No error return*/
    396415        return noerr;
    397         /*}}}*/
    398         }
     416
    399417}

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

    r2771 r2780  
    88struct BamgArgs{
    99
    10         char* geomfile;
     10        int    iso;
     11        double hmin;
     12        double hmax;
     13        double gradation;
     14        double cutoff;
    1115
    1216};

  • issm/trunk/src/m/classes/public/bamg.m

    r2778 r2780  
    1313bamg_mesh=struct();
    1414
    15 %get domainoutline
    16 domainfile=getfieldvalueerr(options,'domain');
    17 resolution=getfieldvalue(options,'resolution',5000);
    18 if ~exist(domainfile,'file')
    19         error(['bamg error message: file ' domainfile ' not found ']);
    20 else
    21         domain=expread(domainfile);
     15% Bamg Geometry parameters {{{1
     16bamg_geometry.NumVertices=0;
     17bamg_geometry.Vertices=zeros(0,3);
     18bamg_geometry.NumEdges=0;
     19bamg_geometry.Edges=zeros(0,3);
     20bamg_geometry.hVertices=zeros(0,1);
     21bamg_geometry.NumSubDomain=0;
     22bamg_geometry.SubDomain=zeros(0,3);
     23if exist(options,'domain'),
     24        domainfile=getfieldvalueerr(options,'domain');
     25        if ~exist(domainfile,'file')
     26                error(['bamg error message: file ' domainfile ' not found ']);
     27        else
     28                domain=expread(domainfile);
    2229
    23         %build geometry
    24         count=0;
    25         bamg_geometry.Vertices=zeros(0,3);
    26         bamg_geometry.Edges=zeros(0,3);
    27         bamg_geometry.hVertices=zeros(0,1);
    28         bamg_geometry.SubDomain=zeros(0,3);
    29         for i=1:length(domain),
    30                 nods=domain(i).nods-1; %the domain are closed 1=end;
    31                 bamg_geometry.Vertices=[bamg_geometry.Vertices; [domain(i).x(1:nods) domain(i).y(1:nods) ones(nods,1)]];
    32                 bamg_geometry.Edges=[bamg_geometry.Edges; [transp(count+1:count+nods) transp([count+2:count+nods count+1])  ones(nods,1)]];
    33                 bamg_geometry.hVertices=[bamg_geometry.hVertices; resolution*ones(nods,1)];
    34                 if i>1,
    35                         bamg_geometry.SubDomain=[2 count+1 1 1];
     30                %build geometry
     31                count=0;
     32
     33                for i=1:length(domain),
     34                        nods=domain(i).nods-1; %the domain are closed 1=end;
     35                        bamg_geometry.Vertices=[bamg_geometry.Vertices; [domain(i).x(1:nods) domain(i).y(1:nods) ones(nods,1)]];
     36                        bamg_geometry.Edges=[bamg_geometry.Edges; [transp(count+1:count+nods) transp([count+2:count+nods count+1])  ones(nods,1)]];
     37                        bamg_geometry.hVertices=[];%[bamg_geometry.hVertices; resolution*ones(nods,1)];
     38                        if i>1,
     39                                clockwise=-1;
     40                                bamg_geometry.SubDomain=[2 count+1 clockwise 1];
     41                        end
     42                        count=count+nods;
    3643                end
    37                 count=count+nods;
    38         end
    39         bamg_geometry.NumVertices=size(bamg_geometry.Vertices,1);
    40         bamg_geometry.NumEdges=size(bamg_geometry.Edges,1);
    41         bamg_geometry.NumSubDomain=size(bamg_geometry.SubDomain,1);
     44                bamg_geometry.NumVertices=size(bamg_geometry.Vertices,1);
     45                bamg_geometry.NumEdges=size(bamg_geometry.Edges,1);
     46                bamg_geometry.NumSubDomain=size(bamg_geometry.SubDomain,1);
     47        end
    4248end
     49%}}}
     50
     51% Bamg Mesh parameters {{{1
     52bamg_mesh.numberofelements=0;
     53bamg_mesh.numberofnodes=0;
     54bamg_mesh.x=[];
     55bamg_mesh.y=[];
     56bamg_mesh.index=zeros(0,3);
     57if (~exist(options,'domain') & md.numberofgrids~=0 & strcmpi(md.type,'2d')),
     58        bamg_mesh.numberofelements=md.numberofelements;
     59        bamg_mesh.numberofnodes=md.numberofgrids;
     60        bamg_mesh.x=md.x;
     61        bamg_mesh.y=md.y;
     62        bamg_mesh.index=md.elements;
     63end
     64%}}}
     65
     66% Bamg Options {{{1
     67bamg_options.iso=getfieldvalue(options,'iso',0);
     68bamg_options.hmin=getfieldvalue(options,'hmin',10^-100);
     69bamg_options.hmax=getfieldvalue(options,'hmax',10^100);
     70bamg_options.gradation=getfieldvalue(options,'gradation',1.2);
     71bamg_options.cutoff=getfieldvalue(options,'cutoff',10^-5);
     72%}}}
    4373
    4474%call Bamg
    45 [elements,x,y]=Bamg(bamg_geometry,bamg_options);
     75[elements,x,y]=Bamg(bamg_mesh,bamg_geometry,bamg_options);
    4676
    4777% plug results onto model

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

    r2778 r2780  
    1313        BamgGeom bamggeom;
    1414
    15         /*inputs: */
    16         char*   geomfile=NULL;
     15        /*Mesh inputs*/
     16        int     numberofnodes;
     17        int     numberofelements;
     18        double* x;
     19        double* y;
     20        double* index;
     21
     22        /*Geom inputs: */
    1723        int     NumVertices;
    1824        double* Vertices=NULL;
    … …  
    2329        double* SubDomain=NULL;
    2430
     31        /*Options inputs*/
     32        int    iso;
     33        double hmin,hmax;
     34        double gradation;
     35        double cutoff;
    2536
    2637        /*Boot module: */
    … …  
    3041        CheckNumMatlabArguments(nlhs,NLHS,nrhs,NRHS,__FUNCT__,&BamgUsage);
    3142
    32         /*process inputs*/
    33         //FetchData(&geomfile,mxGetField(BAMGOPTIONS,0,"fgeom"));
    34 
     43        /*create bamg geometry input*/
    3544        FetchData(&NumVertices,mxGetField(BAMGGEOMETRY,0,"NumVertices"));
     45        bamggeom.NumVertices=NumVertices;
    3646        FetchData(&Vertices,NULL,NULL,mxGetField(BAMGGEOMETRY,0,"Vertices"));
     47        bamggeom.Vertices=Vertices;
    3748        FetchData(&NumEdges,mxGetField(BAMGGEOMETRY,0,"NumEdges"));
     49        bamggeom.NumEdges=NumEdges;
    3850        FetchData(&Edges,NULL,NULL,mxGetField(BAMGGEOMETRY,0,"Edges"));
     51        bamggeom.Edges=Edges;
    3952        FetchData(&hVertices,NULL,NULL,mxGetField(BAMGGEOMETRY,0,"hVertices"));
    40         FetchData(&NumSubDomain,mxGetField(BAMGGEOMETRY,0,"NumSubDomain"));
    41         FetchData(&SubDomain,NULL,NULL,mxGetField(BAMGGEOMETRY,0,"SubDomain"));
    42 
    43         /*create bamg geometry input*/
    44         bamggeom.NumVertices=NumVertices;
    45         bamggeom.Vertices=Vertices;
    46         bamggeom.NumEdges=NumEdges;
    47         bamggeom.Edges=Edges;
    4853        bamggeom.hVertices=hVertices;
    4954        bamggeom.Edges=Edges;
    … …  
    5964        bamggeom.NumRequiredEdges=0;
    6065        bamggeom.RequiredEdges=NULL;
     66        FetchData(&NumSubDomain,mxGetField(BAMGGEOMETRY,0,"NumSubDomain"));
    6167        bamggeom.NumSubDomain=NumSubDomain;
     68        FetchData(&SubDomain,NULL,NULL,mxGetField(BAMGGEOMETRY,0,"SubDomain"));
    6269        bamggeom.SubDomain=SubDomain;
    6370
    6471        /*create bamg mesh input*/
    65         //bamgargs.geomfile=geomfile;
     72        FetchData(&numberofnodes,mxGetField(BAMGMESH,0,"numberofnodes"));
     73        bamgmesh.numberofnodes=numberofnodes;
     74        FetchData(&numberofelements,mxGetField(BAMGMESH,0,"numberofelements"));
     75        bamgmesh.numberofelements=numberofelements;
     76        FetchData(&x,NULL,NULL,mxGetField(BAMGMESH,0,"x"));
     77        bamgmesh.x=x;
     78        FetchData(&y,NULL,NULL,mxGetField(BAMGMESH,0,"y"));
     79        bamgmesh.y=y;
     80        FetchData(&index,NULL,NULL,mxGetField(BAMGMESH,0,"index"));
     81        bamgmesh.index=index;
    6682
    67         /*create bamg mesh input*/
     83        /*create bamg options input*/
     84        FetchData(&iso,mxGetField(BAMGOPTIONS,0,"iso"));
     85        bamgargs.iso=iso;
     86        FetchData(&hmin,mxGetField(BAMGOPTIONS,0,"hmin"));
     87        bamgargs.hmin=hmin;
     88        FetchData(&hmax,mxGetField(BAMGOPTIONS,0,"hmax"));
     89        bamgargs.hmax=hmax;
     90        FetchData(&gradation,mxGetField(BAMGOPTIONS,0,"gradation"));
     91        bamgargs.gradation=gradation;
     92        FetchData(&cutoff,mxGetField(BAMGOPTIONS,0,"cutoff"));
     93        bamgargs.cutoff=cutoff;
    6894
    6995        /*!Generate internal degree of freedom numbers: */
    … …  
    76102
    77103        /*Free ressources: */
    78         //xfree((void**)&geomfile);
    79104        xfree((void**)&Vertices);
    80105        xfree((void**)&Edges);
    … …  
    89114{
    90115        _printf_("\n");
    91         _printf_("   usage: [elements,x,y]=%s(elements,x,y,metric,gradation,splitpbedge,nbv);\n",__FUNCT__);
     116        _printf_("   usage: [elements,x,y]=%s(bamgmesh,bamggeom,bamgoptions,nbv);\n",__FUNCT__);
    92117        _printf_("\n");
    93118}

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

    r2771 r2780  
    1818   
    1919/* serial input macros: */
    20 #define BAMGGEOMETRY (mxArray*)prhs[0]
    21 #define BAMGOPTIONS  (mxArray*)prhs[1]
     20#define BAMGMESH     (mxArray*)prhs[0]
     21#define BAMGGEOMETRY (mxArray*)prhs[1]
     22#define BAMGOPTIONS  (mxArray*)prhs[2]
    2223
    2324/* serial output macros: */
    … …  
    3031#define NLHS  3
    3132#undef NRHS
    32 #define NRHS  2
     33#define NRHS  3
    3334
    3435#endif  /* _BAMG_H */
Note: See TracChangeset for help on using the changeset viewer.