AFEPack
EasyMesh.templates.h
浏览该文件的文档。
00001 
00011 #include "EasyMesh.h"
00012 
00013 template <int DOW>
00014 void TriangleMesh<DOW>::readData(const std::string& filename)
00015 {
00016   int i, j;
00017   int n_node, n_side, n_element;
00018   char text[64];
00019 
00020   std::cout << "Reading easymesh data file ..." << std::endl;   
00021   std::cout << "\treading node data ..." << std::flush;
00022   std::ifstream is((filename + ".n").c_str());
00023   is >> n_node >> n_element >> n_side;
00024   is.getline(text, 64);
00025   this->Mesh<2,DOW>::point().resize(n_node);
00026   this->Mesh<2,DOW>::geometry(0).resize(n_node);
00027   for (i = 0;i < n_node;i ++) {
00028     is >> j
00029        >> this->Mesh<2,DOW>::point(i)
00030        >> this->Mesh<2,DOW>::boundaryMark(0,i);
00031     this->Mesh<2,DOW>::geometry(0,i).index() = j;
00032     this->Mesh<2,DOW>::geometry(0,i).vertex().resize(1,j);
00033     this->Mesh<2,DOW>::geometry(0,i).boundary().resize(1,j);
00034   }
00035   is.close();
00036   std::cout << " OK!" << std::endl;
00037         
00038   std::cout << "\treading side data ..." << std::flush;
00039   is.open((filename + ".s").c_str());
00040   is >> i;
00041   Assert(i == n_side, ExcMeshData("in side file: side number error"));
00042   this->Mesh<2,DOW>::geometry(1).resize(n_side);
00043   for (i = 0;i < n_side;i ++) {
00044     Geometry& g = this->Mesh<2,DOW>::geometry(1,i);
00045     g.vertex().resize(2);
00046     is >> g.index()
00047        >> g.vertex(0) >> g.vertex(1)
00048        >> j >> j
00049        >> this->Mesh<2,DOW>::boundaryMark(1,i);
00050     this->Mesh<2,DOW>::geometry(1,i).boundary() = this->Mesh<2,DOW>::geometry(1,i).vertex();
00051   }
00052   is.close();
00053   std::cout << " OK!" << std::endl;
00054         
00055   std::cout << "\treading element data ..." << std::flush;
00056   is.open((filename + ".e").c_str());
00057   is >> i;
00058   Assert(i == n_element, ExcMeshData("in element file: element number error"));
00059   is >> i;
00060   Assert(i == n_node, ExcMeshData("in element file: node number error"));
00061   is >> i;
00062   Assert(i == n_side, ExcMeshData("in element file: side number error"));
00063   is.getline(text, 64);
00064   this->Mesh<2,DOW>::geometry(2).resize(n_element);
00065   for (i = 0;i < n_element;i ++) {
00066     Geometry& g = this->Mesh<2,DOW>::geometry(2,i);
00067     g.vertex().resize(3);
00068     g.boundary().resize(3);
00069     is >> g.index()
00070        >> g.vertex(0) >> g.vertex(1) >> g.vertex(2)
00071        >> j >> j >> j
00072        >> g.boundary(0) >> g.boundary(1) >> g.boundary(2);
00073     this->Mesh<2,DOW>::boundaryMark(2,i) = 0;
00074 #if (DOW==2)
00075     const Point<DOW>& p0 = this->Mesh<2,DOW>::point(g.vertex(0));
00076     const Point<DOW>& p1 = this->Mesh<2,DOW>::point(g.vertex(1));
00077     const Point<DOW>& p2 = this->Mesh<2,DOW>::point(g.vertex(2));
00078     if ((p1[0] - p0[0])*(p2[1] - p0[1]) - (p1[1] - p0[1])*(p2[0] - p0[0]) < 0) {
00079       std::cerr << "vertices of element " << i 
00080                 << " is not correctly ordered." << std::endl;
00081       j = g.vertex(2);
00082       g.vertex(2) = g.vertex(1);
00083       g.vertex(1) = j;
00084       j = g.boundary(2);
00085       g.boundary(2) = g.boundary(1);
00086       g.boundary(1) = j;
00087     }
00088 #endif
00089   }
00090   is.close();
00091   std::cout << " OK!" << std::endl;
00092 }
00093 
00094 
00095 template <int DOW>
00096 void TriangleMesh<DOW>::writeData(const std::string& filename) const
00097 {
00098   int i, j, k;
00099         
00100   std::cout << "Writing easymesh data file ..." << std::endl;   
00101   std::cout << "\tpreparing data ..." << std::flush;
00102   int n_node = this->Mesh<2,DOW>::n_geometry(0);
00103   int n_side = this->Mesh<2,DOW>::n_geometry(1);
00104   int n_element = this->Mesh<2,DOW>::n_geometry(2);
00105   std::vector<std::vector<int> > side_neighbour(2,
00106                                                 std::vector<int>(n_side, -1));
00107   std::vector<std::vector<int> > element_neighbour(3,
00108                                                    std::vector<int>(n_element, -1));
00109 
00110   for (i = 0;i < n_element;i ++) {
00111     const GeometryBM& the_element = this->Mesh<2,DOW>::geometry(2, i);
00112     for (j = 0;j < 3;j ++) {
00113       k = the_element.boundary(j);
00114       const GeometryBM& the_side = this->Mesh<2,DOW>::geometry(1, k);
00115       if (the_side.vertex(0) == the_element.vertex((j+1)%3)) {
00116         Assert(the_side.vertex(1) == the_element.vertex((j+2)%3), ExcInternalError());
00117         side_neighbour[0][k] = i;
00118       }
00119       else if (the_side.vertex(0) == the_element.vertex((j+2)%3)) {
00120         Assert(the_side.vertex(1) == the_element.vertex((j+1)%3), ExcInternalError());
00121         side_neighbour[1][k] = i;
00122       }
00123     }
00124   }
00125   for (i = 0;i < n_element;i ++) {
00126     const GeometryBM& the_element = this->Mesh<2,DOW>::geometry(2, i);
00127     for (j = 0;j < 3;j ++) {
00128       k = the_element.boundary(j);
00129       const GeometryBM& the_side = this->Mesh<2,DOW>::geometry(1, k);
00130       if (the_side.vertex(0) == the_element.vertex((j+1)%3)) {
00131         element_neighbour[j][i] = side_neighbour[1][k];
00132       }
00133       else if (the_side.vertex(0) == the_element.vertex((j+2)%3)) {
00134         element_neighbour[j][i] = side_neighbour[0][k];
00135       }
00136       else {
00137         Assert(false, ExcInternalError()); // something must be wrong!
00138       }
00139     }
00140   }
00141   std::cout << " OK!" << std::endl;
00142 
00143   std::cout << "\twriting node data ..." << std::flush;
00144   std::ofstream os((filename+".n").c_str());
00145   os.precision(12);
00146   os.setf(std::ios::scientific, std::ios::floatfield);
00147   os << n_node << "\t"
00148      << n_element << "\t"
00149      << n_side << " **(Nnd, Nee, Nsd)**\n";
00150   for (i = 0;i < n_node;i ++) {
00151     j = this->Mesh<2,DOW>::geometry(0, i).vertex(0);
00152     os << i << "\t"
00153        << this->Mesh<2,DOW>::point(j) << "\t"
00154        << this->Mesh<2,DOW>::boundaryMark(0, i) << "\n";
00155   }
00156   os << "---------------------------------------------------------------\n"
00157      << "   n:  x                          y                            \n";
00158   os.close();
00159   std::cout << " OK!" << std::endl;
00160         
00161   std::cout << "\twriting side data ..." << std::flush;
00162   os.open((filename+".s").c_str());
00163   os << n_side << "\n";
00164   for (i = 0;i < n_side;i ++) {
00165     const GeometryBM& the_side = this->Mesh<2,DOW>::geometry(1, i);
00166     os << i << "\t"
00167        << the_side.vertex(0) << "\t"
00168        << the_side.vertex(1) << "\t"
00169        << side_neighbour[0][i] << "\t"
00170        << side_neighbour[1][i] << "\t"
00171        << this->Mesh<2,DOW>::boundaryMark(1, i) << "\n";
00172   }
00173   os << "---------------------------------------------------------------\n"
00174      << "   s:    c      d       ea       eb                            \n";
00175   os.close();
00176   std::cout << " OK!" << std::endl;
00177 
00178   std::cout << "\twriting element data ..." << std::flush;
00179   os.open((filename+".e").c_str());
00180   os << n_element << "\t"
00181      << n_node << "\t"
00182      << n_side << " **(Nee, Nnd, Nsd)**\n";
00183   for (i = 0;i < n_element;i ++) {
00184     const GeometryBM& the_element = this->Mesh<2,DOW>::geometry(2, i);
00185     os << i << "\t"
00186        << the_element.vertex(0) << "\t"
00187        << the_element.vertex(1) << "\t"
00188        << the_element.vertex(2) << "\t"
00189        << element_neighbour[0][i] << "\t"
00190        << element_neighbour[1][i] << "\t"
00191        << element_neighbour[2][i] << "\t"
00192        << the_element.boundary(0) << "\t"
00193        << the_element.boundary(1) << "\t"
00194        << the_element.boundary(2) << "\n";
00195   }
00196   os << "---------------------------------------------------------------\n"
00197      << "   e:  i,   j,   k,   ei,   ej,   ek,   si,   sj,   sk         \n";
00198   os.close();
00199   std::cout << " OK!" << std::endl;
00200 }
00201