Changeset 2772


Ignore:
Timestamp:
01/06/10 14:08:38 (15 years ago)
Author:
Mathieu Morlighem
Message:

now Bamg does not use files anymore (youhouss)

Location:
issm/trunk/src
Files:
6 edited

Legend:

Unmodified
Added
Removed

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

    r2771 r2772  
    8282        double* elements=NULL;
    8383
    84         printf("NumVertices = %i\nNumEdges = %i\n",bamggeom->NumVertices,bamggeom->NumEdges);
    85 
    8684        /*Recover options from inputs: */
    8785        if(splitpbedge)SplitEdgeWith2Boundary=1;
    … …  
    111109
    112110        if(1){
    113                 if (verbosity)
    114                  printf("Construction of a mesh from a given geometry\n");
    115 
    116                 Geometry Gh(bamgargs->geomfile);
     111                if (verbosity>0) printf("Construction of a mesh from a given geometry\n");
     112                Geometry Gh(bamggeom);
    117113                hmin = Max(hmin,Gh.MinimalHmin());
    118114                hmax = Min(hmax,Gh.MaximalHmax());
    … …  
    131127                Triangles Th(nbvx,Gh);
    132128
    133                 /*Build output{{{2*/
     129                /*Build output {{{1*/
    134130                nelout=Th.nbt-Th.NbOutT; //number of triangles - number of external triangles
    135131                nodsout=Th.nbv;
    … …  
    163159
    164160                return noerr;
    165                 /*}}}2*/
     161                /*}}}1*/
     162
    166163        }
    167164

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

    r2740 r2772  
    2828#ifndef _MESH2_H_
    2929#define _MESH2_H_
     30
     31#include "../objects/objects.h"
    3032
    3133#include <stdlib.h>
    … …  
    982984  Real8 MaximalHmax() {return Max(pmax.x-pmin.x,pmax.y-pmin.y);}
    983985  void ReadGeometry(const char * ) ;
     986  void ReadGeometry(BamgGeom* bamggeom);
    984987  void ReadGeometry(MeshIstream & ,const char *)  ;
    985988
    … …  
    987990  Geometry() {EmptyGeometry();}// empty Geometry
    988991  void AfterRead();
     992  Geometry(BamgGeom* bamggeom) {EmptyGeometry();OnDisk=1;ReadGeometry(bamggeom);AfterRead();}
    989993  Geometry(const char * filename) {EmptyGeometry();OnDisk=1;ReadGeometry(filename);AfterRead();}
    990994

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

    r2771 r2772  
    976976
    977977  ReadFromMatlabMesh(elements,x,y,nel,nods,NBV,cutoffradian);
    978 printf("ok2\n");
     978  printf("ok2\n");
    979979
    980980  SetIntCoor();
    … …  
    986986
    987987}
     988
    988989void Geometry::ReadGeometry(const char * filename)
    989990{
    … …  
    995996  MeshIstream f_in (filename);
    996997  ReadGeometry(f_in,filename);
     998}
     999
     1000void Geometry::ReadGeometry(BamgGeom* bamggeom){
     1001
     1002        int verbose=1;
     1003        assert(empty());
     1004        nbiv=nbv=nbvx=0;
     1005        nbe=nbt=nbtx=0;
     1006        NbOfCurves=0;
     1007
     1008        Real8 Hmin = HUGE_VAL;// the infinie value
     1009        Int4 hvertices=0;
     1010        int i,j,k,n;
     1011
     1012        //initialize some variables
     1013        int Version=1,dim=2;
     1014        nbv=bamggeom->NumVertices;
     1015        nbe=bamggeom->NumEdges;
     1016        nbvx = nbv;
     1017        nbiv = nbv;
     1018
     1019        //some checks
     1020        if(verbose>3) printf("Enter ReadGeometry\n");
     1021        if (bamggeom->NumVertices<=0 || bamggeom->Vertices==NULL){
     1022                throw ErrorException(__FUNCT__,exprintf("the domain provided does not contain any vertex"));
     1023        }
     1024        if (bamggeom->NumEdges<=0 || bamggeom->Edges==NULL){
     1025                throw ErrorException(__FUNCT__,exprintf("the domain provided does not contain any edge"));
     1026        }
     1027
     1028        //Vertices
     1029        if (bamggeom->Vertices){
     1030                if(verbose>3) printf("Reading Vertices\n");
     1031                vertices = new GeometricalVertex[nbvx];
     1032                for (i=0;i<nbv;i++) {
     1033                        vertices[i].r.x=(double)bamggeom->Vertices[i*3+0];
     1034                        vertices[i].r.y=(double)bamggeom->Vertices[i*3+1];
     1035                        vertices[i].ReferenceNumber=(Int4)bamggeom->Vertices[i*3+2];
     1036                        vertices[i].DirOfSearch=NoDirOfSearch;
     1037                        vertices[i].color =0;
     1038                        vertices[i].Set();
     1039                }
     1040                //find domain extrema for pmin,pmax
     1041                pmin =  vertices[0].r;
     1042                pmax =  vertices[0].r;
     1043                for (i=0;i<nbv;i++) {
     1044                        pmin.x = Min(pmin.x,vertices[i].r.x);
     1045                        pmin.y = Min(pmin.y,vertices[i].r.y);
     1046                        pmax.x = Max(pmax.x,vertices[i].r.x);
     1047                        pmax.y = Max(pmax.y,vertices[i].r.y);
     1048                }
     1049                R2 DD05 = (pmax-pmin)*0.05;
     1050                pmin -=  DD05;
     1051                pmax +=  DD05;
     1052                coefIcoor= (MaxICoor)/(Max(pmax.x-pmin.x,pmax.y-pmin.y));
     1053                if(coefIcoor <=0){
     1054                        throw ErrorException(__FUNCT__,exprintf("coefIcoor should be positive"));
     1055                }
     1056        }
     1057        else{
     1058                throw ErrorException(__FUNCT__,exprintf("No Vertex provided"));
     1059        }
     1060
     1061        //Edges
     1062        if (bamggeom->Edges){
     1063                int i1,i2;
     1064                R2 zero2(0,0);
     1065                Real4 *len =0;
     1066
     1067                if(verbose>3) printf("Reading Edges\n");
     1068                edges = new GeometricalEdge[nbe];
     1069
     1070                if (!hvertices) {
     1071                        len = new Real4[nbv];
     1072                        for(i=0;i<nbv;i++)
     1073                         len[i]=0;
     1074                }
     1075
     1076                for (i=0;i<nbe;i++){
     1077                        i1=(int)bamggeom->Edges[i*3+0]-1; //-1 for C indexing
     1078                        i2=(int)bamggeom->Edges[i*3+1]-1; //-1 for C indexing
     1079                        edges[i].ref=(Int4)bamggeom->Edges[i*3+2];
     1080                        edges[i].v[0]=  vertices + i1;
     1081                        edges[i].v[1]=  vertices + i2;
     1082                        R2 x12 = vertices[i2].r-vertices[i1].r;
     1083                        Real8 l12=sqrt((x12,x12));
     1084                        edges[i].tg[0]=zero2;
     1085                        edges[i].tg[1]=zero2;
     1086                        edges[i].SensAdj[0] = edges[i].SensAdj[1] = -1;
     1087                        edges[i].Adj[0] = edges[i].Adj[1] = 0;
     1088                        edges[i].flag = 0;
     1089                        if (!hvertices) {
     1090                                vertices[i1].color++;
     1091                                vertices[i2].color++;
     1092                                len[i1] += l12;
     1093                                len[i2] += l12;
     1094                        }
     1095
     1096                        Hmin = Min(Hmin,l12);
     1097                }
     1098
     1099                // definition  the default of the given mesh size
     1100                if (!hvertices){
     1101                        for (i=0;i<nbv;i++)
     1102                         if (vertices[i].color > 0)
     1103                          vertices[i].m=  Metric(len[i] /(Real4) vertices[i].color);
     1104                         else
     1105                          vertices[i].m=  Metric(Hmin);
     1106                        delete [] len;
     1107                }
     1108        }
     1109        else{
     1110                throw ErrorException(__FUNCT__,exprintf("No edges provided"));
     1111        }
     1112
     1113        //hVertices
     1114        if(bamggeom->hVertices){
     1115                if(verbose>3) printf("Reading hVertices\n");
     1116                for (i=0;i< nbv;i++){
     1117                        vertices[i].m=Metric((Real4)bamggeom->hVertices[i]);
     1118                }
     1119        }
     1120        else{
     1121                if(verbose>3) printf("No hVertices found\n");
     1122        }
     1123
     1124        //MetricVertices
     1125        if(bamggeom->MetricVertices){
     1126                if(verbose>3) printf("Reading MetricVertices\n");
     1127                hvertices=1;
     1128                for (i=0;i< nbv;i++) {
     1129                        vertices[i].m = Metric((Real4)bamggeom->MetricVertices[i*3+0],(Real4)bamggeom->MetricVertices[i*3+1],(Real4)bamggeom->MetricVertices[i*3+2]);
     1130                }
     1131        }
     1132        else{
     1133                if(verbose>3) printf("No MetricVertices found\n");
     1134        }
     1135
     1136        //h1h2VpVertices
     1137        if(bamggeom->h1h2VpVertices){
     1138                if(verbose>3) printf("Reading h1h2VpVertices\n");
     1139                Real4 h1,h2,v1,v2;
     1140                hvertices =1;
     1141                for (i=0;i< nbv;i++) {
     1142                        h1=(Real4)bamggeom->MetricVertices[i*4+0];
     1143                        h2=(Real4)bamggeom->MetricVertices[i*4+1];
     1144                        v1=(Real4)bamggeom->MetricVertices[i*4+2];
     1145                        v2=(Real4)bamggeom->MetricVertices[i*4+3];
     1146                        vertices[i].m = Metric(MatVVP2x2(1/(h1*h1),1/(h2*h2),D2(v1,v2)));
     1147                }
     1148        }
     1149        else{
     1150                if(verbose>3) printf("No h1h2VpVertices found\n");
     1151        }
     1152
     1153        //MaximalAngleOfCorner
     1154        if (bamggeom->MaximalAngleOfCorner){
     1155                if(verbose>3) printf("Reading MaximalAngleOfCorner\n");
     1156                MaximalAngleOfCorner=bamggeom->MaximalAngleOfCorner*Pi/180;
     1157        }
     1158        else{
     1159                if(verbose>3) printf("No MaximalAngleOfCorner found\n");
     1160        }
     1161
     1162        //TangentAtEdges
     1163        if (bamggeom->TangentAtEdges){
     1164                if(verbose>3) printf("Reading TangentAtEdges");
     1165                int n,i,j,k;
     1166                R2 tg;
     1167
     1168                n=bamggeom->NumTangentAtEdges;
     1169                for (k=0;k<n;k++) {
     1170                        i=(int)bamggeom->TangentAtEdges[k*4+0]-1; //for C indexing
     1171                        j=(int)bamggeom->TangentAtEdges[k*4+1]-1; //for C indexing
     1172                        tg.x=bamggeom->TangentAtEdges[k*4+2];
     1173                        tg.y=bamggeom->TangentAtEdges[k*4+3];
     1174                        if (j!=0 && j!=1){
     1175                                throw ErrorException(__FUNCT__,exprintf("TangentAtEdges second index should be 1 or 2 only"));
     1176                        }
     1177                        edges[i].tg[j] = tg;
     1178                }
     1179        }
     1180        else{
     1181                if(verbose>3) printf("No TangentAtEdges found\n");
     1182        }
     1183
     1184        //Corners
     1185        if(bamggeom->Corners){
     1186                if(verbose>3) printf("Reading Corners");
     1187                n=bamggeom->NumCorners;
     1188                for (i=0;i<n;i++) {     
     1189                        j=(int)bamggeom->Corners[i]-1; //for C indexing
     1190                        if (j>nbv-1 || j<0){
     1191                                throw ErrorException(__FUNCT__,exprintf("Bad corner definition: should in [0 %i]",nbv));
     1192                        }
     1193                        vertices[j].SetCorner();
     1194                        vertices[j].SetRequired();  }
     1195        }
     1196        else{
     1197                if(verbose>3) printf("No Corners found\n");
     1198        }
     1199
     1200        //RequiredVertices
     1201        if(bamggeom->RequiredVertices){
     1202                if(verbose>3) printf("Reading RequiredVertices");
     1203                n=bamggeom->NumRequiredVertices;
     1204                for (i=0;i<n;i++) {     
     1205                        j=(int)bamggeom->RequiredVertices[i]-1; //for C indexing
     1206                        if (j>nbv-1 || j<0){
     1207                                throw ErrorException(__FUNCT__,exprintf("Bad RequiredVerticess  definition: should in [0 %i]",nbv));
     1208                        }
     1209                        vertices[j].SetRequired();  }
     1210        }
     1211        else{
     1212                if(verbose>3) printf("No RequiredVertices found\n");
     1213        }
     1214
     1215        //RequiredEdges
     1216        if(bamggeom->RequiredEdges){
     1217                if(verbose>3) printf("Reading RequiredEdges");
     1218                n=bamggeom->NumRequiredEdges;
     1219                for (i=0;i<n;i++) {     
     1220                        j=(int)bamggeom->RequiredEdges[i]-1; //for C indexing
     1221                        if (j>nbe-1 || j<0){
     1222                                throw ErrorException(__FUNCT__,exprintf("Bad RequiredEdges definition: should in [0 %i]",nbv));
     1223                        }
     1224                        edges[j].SetRequired();  }
     1225        }
     1226        else{
     1227                if(verbose>3) printf("No RequiredEdges found\n");
     1228        }
     1229
     1230        //SubDomain
     1231        if(bamggeom->SubDomain){
     1232                Int4 i0,i1,i2,i3;
     1233                if(verbose>3) printf("Reading SubDomain");
     1234                NbSubDomains=bamggeom->NumSubDomain;
     1235                subdomains = new GeometricalSubDomain[NbSubDomains];
     1236                for (i=0;i<NbSubDomains;i++) {
     1237                        i0=(int)bamggeom->SubDomain[i*4+0];
     1238                        i1=(int)bamggeom->SubDomain[i*4+1];
     1239                        i2=(int)bamggeom->SubDomain[i*4+2];
     1240                        i3=(int)bamggeom->SubDomain[i*4+3];
     1241                        if (i0!=2) throw ErrorException(__FUNCT__,exprintf("Bad Subdomain definition: first number should be 2"));
     1242                        if (i1>nbe || i1<=0) throw ErrorException(__FUNCT__,exprintf("Bad Subdomain definition: second number should in [1 %i]",nbe));
     1243                        subdomains[i].edge=edges + (i1-1);
     1244                        subdomains[i].sens = (int) i2;
     1245                        subdomains[i].ref = i3;
     1246                }
     1247        }
     1248        else{
     1249                if(verbose>3) printf("No SubDomain found\n");
     1250        }
    9971251}
    9981252

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

    r2740 r2772  
    4242SaveMetricInterpole  LastMetricInterpole;
    4343
    44 void ReductionSimultanee( MetricAnIso M1,  MetricAnIso M2,double & l1,double & l2, D2xD2 & V)
    45 {
     44void ReductionSimultanee( MetricAnIso M1,  MetricAnIso M2,double & l1,double & l2, D2xD2 & V) {
    4645  double a11=M1.a11,a21=M1.a21,a22=M1.a22;
    4746  double b11=M2.a11,b21=M2.a21,b22=M2.a22;
    … …  
    108107
    109108MetricAnIso Intersection(const MetricAnIso M1,const MetricAnIso M2) ;
    110 MetricAnIso Intersection(const MetricAnIso M1,const MetricAnIso M2)
    111 {
     109MetricAnIso Intersection(const MetricAnIso M1,const MetricAnIso M2) {
    112110      D2xD2 M;
    113111      double l1,l2;
    … …  
    121119}
    122120
    123 MetricAnIso::MetricAnIso(const Real8  a[3],const  MetricAnIso m0,
    124            const  MetricAnIso m1,const  MetricAnIso m2 )
    125 {
     121MetricAnIso::MetricAnIso(const Real8  a[3],const  MetricAnIso m0, const  MetricAnIso m1,const  MetricAnIso m2 ){
    126122  MetricAnIso mab(a[0]*m0.a11 + a[1]*m1.a11 + a[2]*m2.a11,
    127123                  a[0]*m0.a21 + a[1]*m1.a21 + a[2]*m2.a21,
    … …  
    141137}
    142138
    143  MetricAnIso::MetricAnIso( Real8  a,const  MetricAnIso ma,
    144                            Real8  b,const  MetricAnIso mb)
    145 {
     139 MetricAnIso::MetricAnIso( Real8  a,const  MetricAnIso ma, Real8  b,const  MetricAnIso mb) {
    146140  MetricAnIso mab(a*ma.a11+b*mb.a11,a*ma.a21+b*mb.a21,a*ma.a22+b*mb.a22);
    147141  MatVVP2x2 vab(mab);
    … …  
    160154
    161155
    162  MatVVP2x2::MatVVP2x2(const MetricAnIso M)
    163 {
     156 MatVVP2x2::MatVVP2x2(const MetricAnIso M){
    164157  double a11=M.a11,a21=M.a21,a22=M.a22;
    165158  const double eps = 1.e-5;

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

    r2771 r2772  
    99
    1010        int     NumVertices;
    11         double* Vertices;
     11        double* Vertices;
     12
    1213        int     NumEdges;
    1314        double* Edges;
     15
    1416        double* hVertices;
     17        double* MetricVertices;
     18        double* h1h2VpVertices;
    1519
     20        double  MaximalAngleOfCorner;
     21
     22        int     NumTangentAtEdges;
     23        double* TangentAtEdges;
     24
     25        int     NumCorners;
     26        double* Corners;
     27
     28        int     NumRequiredVertices;
     29        double* RequiredVertices;
     30
     31        int     NumRequiredEdges;
     32        double* RequiredEdges;
     33
     34        int     NumSubDomain;
     35        double* SubDomain;
    1636};
    1737#endif

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

    r2771 r2772  
    3838
    3939        /*create bamg geometry input*/
    40         bamgargs.geomfile=geomfile;
    41 
    4240        bamggeom.NumVertices=NumVertices;
    4341        bamggeom.Vertices=Vertices;
    … …  
    4543        bamggeom.Edges=Edges;
    4644        bamggeom.hVertices=hVertices;
     45        bamggeom.Edges=Edges;
     46        bamggeom.MetricVertices=NULL;
     47        bamggeom.h1h2VpVertices=NULL;
     48        bamggeom.MaximalAngleOfCorner=10;
     49        bamggeom.NumTangentAtEdges=0;
     50        bamggeom.TangentAtEdges=NULL;
     51        bamggeom.NumCorners=0;
     52        bamggeom.Corners=NULL;
     53        bamggeom.NumRequiredVertices=0;
     54        bamggeom.RequiredVertices=NULL;
     55        bamggeom.NumRequiredEdges=0;
     56        bamggeom.RequiredEdges=NULL;
     57        bamggeom.NumSubDomain=0;
     58        bamggeom.SubDomain=NULL;
     59
     60        /*create bamg mesh input*/
     61        bamgargs.geomfile=geomfile;
     62
     63        /*create bamg mesh input*/
    4764
    4865        /*!Generate internal degree of freedom numbers: */
Note: See TracChangeset for help on using the changeset viewer.