Changeset 2797
- Timestamp:
- 01/12/10 11:17:20 (15 years ago)
- Location:
- issm/trunk/src/c/Bamgx
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Bamgx/Mesh2.h
r2796 r2797 870 870 GeometricalEdge * ProjectOnCurve( Edge & AB, Vertex & A, Vertex & B,Real8 theta, 871 871 Vertex & R,VertexOnEdge & BR,VertexOnGeom & GR); 872 873 874 void WriteElements(std::ostream& f,Int4 * reft ,Int4 nbInT) const;875 876 872 877 873 Int4 Number(const Triangle & t) const { return &t - triangles;} … … 907 903 int UnCrack(); 908 904 909 #ifdef DEBUG910 void inline Check();911 #endif912 #ifdef DRAWING913 void Draw() const ;914 void InitDraw() const ;915 void inquire() ;916 #endif917 905 friend std::ostream& operator <<(std::ostream& f, const Triangles & Th); 918 void Write(const char * filename);919 906 void ConsGeometry(Real8 =-1.0,int *equiedges=0); // construct a geometry if no geo 920 907 void FillHoleInMesh() ; … … 924 911 void GeomToTriangles0(Int4 nbvx);// the real constructor mesh generator 925 912 void PreInit(Int4,char * =0 ); 926 //927 void Write_nop5(OFortranUnFormattedFile * f,928 Int4 &lnop5,Int4 &nef,Int4 &lgpdn,Int4 ndsr) const ;929 930 913 931 914 }; // End Class Triangles … … 970 953 Real8 MinimalHmin() {return 2.0/coefIcoor;} 971 954 Real8 MaximalHmax() {return Max(pmax.x-pmin.x,pmax.y-pmin.y);} 972 void ReadGeometry(const char * ) ;973 955 void ReadGeometry(BamgGeom* bamggeom, BamgOpts* bamgopts); 974 void ReadGeometry(MeshIstream & ,const char *) ;975 976 956 void EmptyGeometry(); 977 957 Geometry() {EmptyGeometry();}// empty Geometry 978 958 void AfterRead(); 979 959 Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts) {EmptyGeometry();OnDisk=1;ReadGeometry(bamggeom,bamgopts);AfterRead();} 980 Geometry(const char * filename) {EmptyGeometry();OnDisk=1;ReadGeometry(filename);AfterRead();} 981 982 void ReadMetric(const char *,Real8 hmin,Real8 hmax,Real8 coef); 960 983 961 const GeometricalVertex & operator[] (Int4 i) const { return vertices[i];}; 984 962 GeometricalVertex & operator[](Int4 i) { return vertices[i];}; … … 997 975 GeometricalEdge * Contening(const R2 P, GeometricalEdge * start) const; 998 976 friend std::ostream& operator <<(std::ostream& f, const Geometry & Gh); 999 void Write(const char * filename);1000 977 void WriteGeometry(BamgGeom* bamggeom, BamgOpts* bamgopts); 1001 1002 978 }; // End Class Geometry 1003 979 -
issm/trunk/src/c/Bamgx/MeshRead.cpp
r2795 r2797 529 529 } 530 530 } 531 532 void Geometry::ReadGeometry(MeshIstream & f_in,const char * filename)533 {534 long int verbosity=2;535 if(verbosity>1)536 cout << " -- ReadGeometry " << filename << endl;537 assert(empty());538 nbiv=nbv=nbvx=0;539 nbe=nbt=nbtx=0;540 NbOfCurves=0;541 // BeginOfCurve=0;542 name=new char [strlen(filename)+1];543 strcpy(name,filename);544 Real8 Hmin = HUGE_VAL;// the infinie value545 // Real8 MaximalAngleOfCorner = 10*Pi/180; ;546 Int4 hvertices =0;547 Int4 i;548 Int4 Version,dim=0;549 int field=0;550 int showfield=0;551 int NbErr=0;552 553 while (f_in.cm())554 {555 field=0;556 // warning ruse for on allocate fiedname at each time557 char fieldname[256];558 f_in.cm() >> fieldname ;559 f_in.err();560 if(f_in.eof()) break;561 // cout << fieldname << " line " << LineNumber <<endl;562 if (!strcmp(fieldname,"MeshVersionFormatted") )563 f_in >> Version ;564 else if (!strcmp(fieldname,"End"))565 break;566 else if (!strcmp(fieldname,"end"))567 break;568 else if (!strcmp(fieldname,"Dimension"))569 {570 f_in >> dim ;571 if(verbosity>5)572 cout << " Geom Record Dimension dim = " << dim << endl;573 assert(dim ==2);574 }575 else if (!strcmp(fieldname,"hVertices"))576 {577 if (nbv <=0) {578 cerr<<"Error: the field Vertex is not found before hVertices " << filename<<endl;579 NbErr++;}580 if(verbosity>5)581 cout << " Geom Record hVertices nbv=" << nbv << endl;582 hvertices =1;583 for (i=0;i< nbv;i++)584 {585 Real4 h;586 f_in >> h ;587 vertices[i].m = Metric(h);588 }589 }590 else if (!strcmp(fieldname,"MetricVertices"))591 { hvertices =1;592 if (nbv <=0) {593 cerr<<"Error: the field Vertex is not found before MetricVertices " << filename<<endl;594 NbErr++;}595 if(verbosity>5)596 cout << " Geom Record MetricVertices nbv =" << nbv << endl;597 for (i=0;i< nbv;i++)598 {599 Real4 a11,a21,a22;600 f_in >> a11 >> a21 >> a22 ;601 vertices[i].m = Metric(a11,a21,a22);602 }603 }604 else if (!strcmp(fieldname,"h1h2VpVertices"))605 { hvertices =1;606 if (nbv <=0) {607 cerr<<"Error: the field Vertex is not found before h1h2VpVertices " << filename<<endl;608 NbErr++;}609 if(verbosity>5)610 cout << " Geom Record h1h2VpVertices nbv=" << nbv << endl;611 612 for (i=0;i< nbv;i++)613 {614 Real4 h1,h2,v1,v2;615 f_in >> h1 >> h2 >>v1 >>v2 ;616 vertices[i].m = Metric(MatVVP2x2(1/(h1*h1),1/(h2*h2),D2(v1,v2)));617 }618 }619 else if (!strcmp(fieldname,"Vertices"))620 {621 assert(dim ==2);622 f_in >> nbv ;623 // if(LineError) break;624 nbvx = nbv;625 626 vertices = new GeometricalVertex[nbvx];627 if(verbosity>5)628 cout << " Geom Record Vertices nbv = " << nbv << "vertices = " << vertices<<endl;629 assert(nbvx >= nbv);630 nbiv = nbv;631 for (i=0;i<nbv;i++) {632 f_in >> vertices[i].r.x ;633 // if(LineError) break;634 f_in >> vertices[i].r.y ;635 // if(LineError) break;636 f_in >> vertices[i].ReferenceNumber ;637 vertices[i].DirOfSearch=NoDirOfSearch;638 // vertices[i].m.h = 0;639 vertices[i].color =0;640 vertices[i].Set();}641 // if(LineError) break;642 pmin = vertices[0].r;643 pmax = vertices[0].r;644 // recherche des extrema des vertices pmin,pmax645 for (i=0;i<nbv;i++) {646 pmin.x = Min(pmin.x,vertices[i].r.x);647 pmin.y = Min(pmin.y,vertices[i].r.y);648 pmax.x = Max(pmax.x,vertices[i].r.x);649 pmax.y = Max(pmax.y,vertices[i].r.y);650 }651 652 R2 DD05 = (pmax-pmin)*0.05;653 pmin -= DD05;654 pmax += DD05;655 656 coefIcoor= (MaxICoor)/(Max(pmax.x-pmin.x,pmax.y-pmin.y));657 assert(coefIcoor >0);658 if (verbosity>5) {659 cout << " Geom: min="<< pmin << "max ="<< pmax << " hmin = " << MinimalHmin() << endl;}660 }661 else if(!strcmp(fieldname,"MaximalAngleOfCorner")||!strcmp(fieldname,"AngleOfCornerBound"))662 {663 f_in >> MaximalAngleOfCorner;664 665 if(verbosity>5)666 cout << " Geom Record MaximalAngleOfCorner " << MaximalAngleOfCorner <<endl;667 MaximalAngleOfCorner *= Pi/180;668 }669 else if (!strcmp(fieldname,"Edges"))670 {671 if (nbv <=0) {672 cerr<<"Error: the field edges is not found before MetricVertices " << filename<<endl;673 NbErr++;}674 else675 {676 int i1,i2;677 R2 zero2(0,0);678 f_in >> nbe ;679 680 edges = new GeometricalEdge[nbe];681 if(verbosity>5)682 cout << " Record Edges: Nb of Edge " << nbe <<endl;683 assert(edges);684 assert (nbv >0);685 Real4 *len =0;686 if (!hvertices)687 {688 len = new Real4[nbv];689 for(i=0;i<nbv;i++)690 len[i]=0;691 }692 693 for (i=0;i<nbe;i++)694 {695 f_in >> i1 >> i2 >> edges[i].ref ;696 697 i1--;i2--; // for C index698 edges[i].v[0]= vertices + i1;699 edges[i].v[1]= vertices + i2;700 R2 x12 = vertices[i2].r-vertices[i1].r;701 Real8 l12=sqrt((x12,x12));702 edges[i].tg[0]=zero2;703 edges[i].tg[1]=zero2;704 edges[i].SensAdj[0] = edges[i].SensAdj[1] = -1;705 edges[i].Adj[0] = edges[i].Adj[1] = 0;706 edges[i].flag = 0;707 if (!hvertices)708 {709 vertices[i1].color++;710 vertices[i2].color++;711 len[i1] += l12;712 len[i2] += l12;713 }714 715 Hmin = Min(Hmin,l12);716 }717 // definition the default of the given mesh size718 if (!hvertices)719 {720 for (i=0;i<nbv;i++)721 if (vertices[i].color > 0)722 vertices[i].m= Metric(len[i] /(Real4) vertices[i].color);723 else724 vertices[i].m= Metric(Hmin);725 delete [] len;726 727 if(verbosity>3)728 cout << " Geom Hmin " << Hmin << endl;729 }730 731 }732 }733 else if (!strcmp(fieldname,"EdgesTangence") ||!strcmp(fieldname,"TangentAtEdges") )734 {735 int n,i,j,k;736 R2 tg;737 f_in >> n ;738 739 if(verbosity>5)740 cout << " Record TangentAtEdges: Nb of Edge " << n <<endl;741 742 for (k=0;k<n;k++)743 {744 f_in >> i >> j ;745 f_in >> tg.x >> tg.y ;746 assert( i <= nbe );747 assert( i > 0 );748 assert ( j == 1 || j==2 );749 i--;j--;// for C index750 edges[i].tg[j] = tg;751 }752 }753 else if (!strcmp(fieldname,"Corners"))754 {755 int i,j,n;756 f_in >> n ;757 if(verbosity>5)758 cout << " Record Corner: Nb of Corner " << n <<endl;759 760 for (i=0;i<n;i++) {761 f_in >> j ;762 assert( j <= nbv );763 assert( j > 0 );764 j--;765 vertices[j].SetCorner();766 vertices[j].SetRequired(); }767 }768 else if (!strcmp(fieldname,"RequiredVertices"))769 {770 int i,j,n;771 f_in >> n ;772 773 for (i=0;i<n;i++) {774 f_in >> j ;775 assert( j <= nbv );776 assert( j > 0 );777 j--;778 vertices[j].SetRequired(); }779 }780 else if (!strcmp(fieldname,"RequiredEdges"))781 {782 int i,j,n;783 f_in >> n ;784 785 for (i=0;i<n;i++) {786 f_in >> j ;787 assert( j <= nbe );788 assert( j > 0 );789 j--;790 edges[j].SetRequired(); }791 }792 else if (!strcmp(fieldname,"SubDomain") || !strcmp(fieldname,"SubDomainFromGeom"))793 {794 f_in >> NbSubDomains ;795 if (NbSubDomains>0)796 {797 subdomains = new GeometricalSubDomain[ NbSubDomains];798 Int4 i0,i1,i2,i3;799 for (i=0;i<NbSubDomains;i++)800 {801 f_in >> i0 >>i1802 >> i2 >>i3 ;803 804 assert(i0 == 2);805 assert(i1<=nbe && i1>0);806 subdomains[i].edge=edges + (i1-1);807 subdomains[i].sens = (int) i2;808 subdomains[i].ref = i3;809 }810 }811 }812 else813 { // unkown field814 field = ++showfield;815 if(showfield==1) // just to show one time816 if (verbosity>3)817 cout << " Warning we skip the field " << fieldname << " at line " << f_in.LineNumber << endl;818 }819 showfield=field; // just to show one time820 } // while !eof()821 // generation de la geometrie822 // 1 construction des aretes823 // construire des aretes en chaque sommets824 825 if (nbv <=0) {826 cerr<<"Error: the field Vertex is not found in " << filename<<endl;827 NbErr++;}828 if(nbe <=0) {829 cerr <<"Error: the field Edges is not found in "<< filename<<endl830 ;NbErr++;}831 if(NbErr) MeshError(1);832 833 834 }835 836 837 531 } // end of namespace bamg -
issm/trunk/src/c/Bamgx/MeshWrite.cpp
r2796 r2797 215 215 } 216 216 217 void Geometry::Write(const char * filename)218 {219 long int verbosity=0;220 221 ofstream f(filename);222 if (f)223 {224 if(verbosity>1)225 cout << " -- write geometry in file " << filename << endl;226 if (name) delete name;227 name = new char[strlen(filename)+1];228 strcpy(name,filename);229 OnDisk =1;230 f << *this;231 }232 }233 234 217 void Geometry::WriteGeometry(BamgGeom* bamggeom, BamgOpts* bamgopts){ 235 218 … … 391 374 } 392 375 393 ostream& operator <<(ostream& f, const Geometry & Gh)394 {395 Int4 NbCorner=0;396 {397 f << "MeshVersionFormatted 0" <<endl;398 f << "\nDimension\n" << 2 << endl;399 // f << "\nIdentifier\n" ;400 // WriteStr(f,Gh.identity);401 // f <<endl;402 }403 int nbreqv=0;404 {405 406 f.precision(12);407 f << "\nVertices\n" << Gh.nbv <<endl;408 for (Int4 i=0;i<Gh.nbv;i++)409 {410 GeometricalVertex & v = Gh.vertices[i];411 if (v.Required()) nbreqv++;412 f << v.r.x << " " << v.r.y << " " << v.ref() << endl;413 if (v.Corner()) NbCorner++;414 }415 }416 417 int nbcracked=0;418 419 {420 int nbreq=0;421 f << "\nEdges\n"<< Gh.nbe << endl;422 for(Int4 ie=0;ie<Gh.nbe;ie++)423 {424 425 GeometricalEdge & e = Gh.edges[ie];426 if (e.Required()) nbreq++;427 if (e.Cracked()) {428 Int4 ie1 = Gh.Number(e.link);429 if (ie <= ie1) ++nbcracked;}430 f << Gh.Number(e[0])+1 << " " << Gh.Number(e[1])+1;431 f << " " << e.ref <<endl;432 }433 434 if (nbcracked)435 {436 f << "\nCrackedEdges\n"<< nbcracked<< endl;437 for(Int4 ie=0;ie<Gh.nbe;ie++)438 {439 GeometricalEdge & e = Gh.edges[ie];440 if (e.Cracked()) {441 Int4 ie1 = Gh.Number(e.link);442 if (ie <= ie1) f << ie+1 << " " << ie1+1<< endl;443 }444 }445 }446 if(nbreq)447 {448 f << "\nRequiredEdges\n"<< nbreq<< endl;449 for(Int4 ie=0;ie<Gh.nbe;ie++)450 {451 GeometricalEdge & e = Gh.edges[ie];452 if (e.Required())453 f << ie+1 << endl;454 }455 }456 457 458 459 }460 461 f << "\nAngleOfCornerBound\n"462 << Gh.MaximalAngleOfCorner*180/Pi << endl;463 if (NbCorner)464 {465 f << "\nCorners\n" << NbCorner << endl;466 for (Int4 i=0,j=0;i<Gh.nbv;i++)467 {468 GeometricalVertex & v = Gh.vertices[i];469 if (v.Corner())470 j++,f << Gh.Number(v)+1 << (j % 5 ? ' ' : '\n');471 }472 473 474 }475 476 if(nbreqv)477 {478 f << "\nRequiredVertices\n"<< nbreqv<< endl;479 for (Int4 j=0,i=0;i<Gh.nbv;i++)480 {481 GeometricalVertex & v = Gh.vertices[i];482 483 if (v.Required())484 j++,f << i+1 << (j % 5 ? ' ' : '\n');485 }486 f << endl;487 }488 489 {490 Int4 i;491 f << "\nSubDomainFromGeom\n" ;492 f << Gh.NbSubDomains<< endl;493 for (i=0;i<Gh.NbSubDomains;i++)494 f << "2 " << Gh.Number(Gh.subdomains[i].edge)+1 << " " << Gh.subdomains[i].sens495 << " " << Gh.subdomains[i].ref << endl;496 }497 {498 Int4 n=0,i;499 500 for(i=0;i< Gh.nbe;i++)501 {502 if(Gh.edges[i].TgA() && Gh.edges[i][0].Corner() )503 n++;504 if(Gh.edges[i].TgB() && Gh.edges[i][1].Corner() )505 n++;506 }507 if (n) {508 f << "TangentAtEdges " << n << endl;509 for(i=0;i< Gh.nbe;i++)510 {511 if (Gh.edges[i].TgA() && Gh.edges[i][0].Corner() )512 f << i+1 << " 1 " << Gh.edges[i].tg[0].x513 << " " << Gh.edges[i].tg[0].y << endl;514 if (Gh.edges[i].TgB() && Gh.edges[i][1].Corner() )515 f << i+1 << " 2 " << Gh.edges[i].tg[1].x516 << " " << Gh.edges[i].tg[1].y << endl;517 }518 519 }}520 // f << " Not Yet Implemented" << endl;521 522 return f;523 }524 525 526 376 } // end of namespace bamg -
issm/trunk/src/c/Bamgx/Metric.cpp
r2790 r2797 1023 1023 } 1024 1024 1025 void Geometry::ReadMetric(const char * fmetrix,Real8 hmin=1.0e-30,Real8 hmax=1.0e30,Real8 coef=1)1026 {1027 long int verbosity=0;1028 1029 hmin = Max(hmin,MinimalHmin());1030 MeshIstream f_metrix(fmetrix);1031 Int4 k,j;1032 f_metrix >> k >> j ;1033 if(verbosity>1)1034 cout << " -- ReadMetric " << fmetrix1035 << ", coef = " << coef1036 << ", hmin = " << hmin1037 << ", hmax = " << hmax1038 << ( (j == 1)? " Iso " : " AnIso " ) << endl;1039 1040 if (k != nbv || !(j == 1 || j == 3)) {1041 cerr << " Error Pb metrix " << k << " <> "1042 << nbv << " or 1 or 3 <> " << j << endl;1043 MeshError(1003);}1044 1045 1046 // Int4 nberr = 0;1047 for (Int4 iv=0;iv<nbv;iv++)1048 {1049 Real8 h;1050 if (j == 1)1051 {1052 f_metrix >> h ;1053 vertices[iv].m=Metric(Max(hmin,Min(hmax, h*coef)));1054 }1055 else if (j==3)1056 {1057 Real8 a,b,c;1058 f_metrix >> a >> b >> c ;1059 MetricAnIso M(a,b,c);1060 MatVVP2x2 Vp(M/coef);1061 Vp.Maxh(hmin);1062 Vp.Minh(hmax);1063 vertices[iv].m = Vp;1064 }1065 }1066 1067 }1068 1069 1070 1025 Real8 LengthInterpole(const MetricAnIso Ma,const MetricAnIso Mb, R2 AB) 1071 1026 {
Note:
See TracChangeset
for help on using the changeset viewer.