AFEPack
|
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