AFEPack
|
00001 00011 template <> 00012 HGeometry<0,DOW>::HGeometry() 00013 : index(0), bmark(0) {} 00014 00015 00016 template <> 00017 HGeometry<1,DOW>::HGeometry() 00018 : index(0), 00019 vertex(n_vertex, NULL), 00020 parent(NULL), 00021 child(n_child, NULL), 00022 bmark(0) {} 00023 00024 00025 00026 template <> 00027 HGeometry<2,DOW>::HGeometry() 00028 : index(0), 00029 vertex(n_vertex, NULL), 00030 boundary(n_boundary, NULL), 00031 parent(NULL), 00032 child(n_child, NULL), 00033 bmark(0) {} 00034 00035 00036 00037 template <> 00038 HGeometry<3,DOW>::HGeometry() 00039 : index(0), 00040 vertex(n_vertex, NULL), 00041 boundary(n_boundary, NULL), 00042 parent(NULL), 00043 child(n_child, NULL), 00044 bmark(0) {} 00045 00046 00047 template <> 00048 bool HGeometry<1,DOW>::isRefined() const 00049 { 00050 return (child.size() > 0 && child[0] != NULL); 00051 } 00052 00053 00054 00055 template <> 00056 bool HGeometry<2,DOW>::isRefined() const 00057 { 00058 return (child.size() > 0 && child[0] != NULL); 00059 } 00060 00061 00062 00063 template <> 00064 bool HGeometry<3,DOW>::isRefined() const 00065 { 00066 return (child.size() > 0 && child[0] != NULL); 00067 } 00068 00069 00070 00071 template <> 00072 void HGeometry<1,DOW>::refine() 00073 { 00074 if (isRefined()) return; 00075 00076 // add the midpoint of the edge at first 00077 HGeometry<0,DOW> * new_point = new HGeometry<0,DOW>(); 00078 Assert(new_point != NULL, ExcOutOfMemory()); 00079 if (mid_point == NULL) { 00080 *dynamic_cast<afepack::Point<DOW> *>(new_point) = midpoint(*vertex[0], *vertex[1]); 00081 } 00082 else { 00083 (*mid_point)(*vertex[0], *vertex[1], bmark, *new_point); 00084 } 00085 new_point->bmark = bmark; 00086 00087 // construct the 0th child for the edge 00088 child[0] = new HGeometry<1,DOW>(); 00089 Assert(child[0] != NULL, ExcOutOfMemory()); 00090 child[0]->parent = this; 00091 child[0]->vertex[0] = vertex[0]; 00092 child[0]->vertex[1] = new_point; 00093 child[0]->bmark = bmark; 00094 00095 // construct the 1st child for the edge 00096 child[1] = new HGeometry<1,DOW>(); 00097 Assert(child[1] != NULL, ExcOutOfMemory()); 00098 child[1]->parent = this; 00099 child[1]->vertex[0] = new_point; 00100 child[1]->vertex[1] = vertex[1]; 00101 child[1]->bmark = bmark; 00102 } 00103 00104 00105 00106 template <> 00107 void HGeometry<2,DOW>::refine() 00108 { 00109 int i, ii[]={0, 1, 2, 0, 1, 2}; 00110 if (isRefined()) return; 00111 00112 // refine its edge at first 00113 for (i = 0;i < 3;i ++) 00114 boundary[i]->refine(); 00115 00116 // add the three new edge in the triangle 00117 HGeometry<0,DOW> * edge_midpoint[3]; 00118 for (i = 0;i < 3;i ++) { 00119 edge_midpoint[i] = boundary[i]->child[0]->vertex[1]; 00120 } 00121 HGeometry<1,DOW> * new_edge[3]; 00122 for (i = 0;i < 3;i ++) { 00123 new_edge[i] = new HGeometry<1,DOW>(); 00124 Assert(new_edge[i] != NULL, ExcOutOfMemory()); 00125 new_edge[i]->vertex[0] = edge_midpoint[ii[i+1]]; 00126 new_edge[i]->vertex[1] = edge_midpoint[ii[i+2]]; 00127 new_edge[i]->bmark = bmark; 00128 } 00129 00130 // construct the 0th to 2nd child for the triangle 00131 for (i = 0;i < 3;i ++) { 00132 if (typeid(*this) == typeid(HGeometry<2, DOW>)) 00133 child[i] = (HGeometry<2,DOW> *) new HGeometry<2, DOW>(); 00134 else 00135 child[i] = new HGeometry<2,DOW>(); 00136 Assert(child[i] != NULL, ExcOutOfMemory()); 00137 child[i]->parent = this; 00138 child[i]->vertex[0] = vertex[i]; 00139 child[i]->vertex[1] = edge_midpoint[ii[i+2]]; 00140 child[i]->vertex[2] = edge_midpoint[ii[i+1]]; 00141 child[i]->boundary[0] = new_edge[i]; 00142 if (boundary[ii[i+1]]->vertex[0] == vertex[i]) 00143 child[i]->boundary[1] = boundary[ii[i+1]]->child[0]; 00144 else 00145 child[i]->boundary[1] = boundary[ii[i+1]]->child[1]; 00146 if (boundary[ii[i+2]]->vertex[0] == vertex[i]) 00147 child[i]->boundary[2] = boundary[ii[i+2]]->child[0]; 00148 else 00149 child[i]->boundary[2] = boundary[ii[i+2]]->child[1]; 00150 child[i]->bmark = bmark; 00151 } 00152 00153 // construct the 3rd child, the center child, for the triangle 00154 child[3] = new HGeometry<2,DOW>(); 00155 Assert(child[3] != NULL, ExcOutOfMemory()); 00156 child[3]->parent = this; 00157 child[3]->vertex[0] = edge_midpoint[0]; 00158 child[3]->vertex[1] = edge_midpoint[1]; 00159 child[3]->vertex[2] = edge_midpoint[2]; 00160 child[3]->boundary[0] = new_edge[0]; 00161 child[3]->boundary[1] = new_edge[1]; 00162 child[3]->boundary[2] = new_edge[2]; 00163 child[3]->bmark = bmark; 00164 } 00165 00166 00167 00168 template <> 00169 void HGeometry<3,DOW>::refine() 00170 { 00171 int i, ii[]={0, 1, 2, 0, 1, 2, 0, 1, 2}; 00172 if (isRefined()) return; 00173 00174 // refine its surfaces at first. After the surfaces of the tetrahedron are refined, 00175 // all those edges of the tetrahedron are refined, too. 00176 for (i = 0;i < HGeometry<3,DOW>::n_boundary;i ++) 00177 boundary[i]->refine(); 00178 00179 // retrieve the information of those midpoints of the edges of the tetrahedron. 00180 // Because the data of the surfaces of the tetrahedron are not unique for the 00181 // tetrahedron, it's not so simple to get the information of those midpoints. 00182 // A lot of codes are written to get the relationship detail between the surfaces 00183 // and the tetrahedron. Be patiently and let's do it step by step, 00184 HGeometry<0,DOW> * edge_midpoint[6]; 00185 HGeometry<2,DOW> * surface_child[4][4]; 00186 00187 { 00188 // for the first surface of the tetrahedron, which is on the bottom of the tetrahedron, 00189 // we should judge the pose of the triangle at first. 00190 // boundary[0]: the first surface 00191 // boundary[0]->vertex[i]: the i-th vertex of the triangle 00192 int start, step; 00193 for (start = 0;start < 3; start ++) 00194 if (boundary[0]->vertex[start] == vertex[1]) 00195 break; 00196 Assert (start < 3,ExcInternalError()); 00197 if (boundary[0]->vertex[ii[start+1]] == vertex[2]) 00198 step = 1; 00199 else { 00200 Assert(boundary[0]->vertex[ii[start+1]] == vertex[3], ExcInternalError()); 00201 Assert(boundary[0]->vertex[ii[start+2]] == vertex[2], ExcInternalError()); 00202 step = 2; 00203 } 00204 edge_midpoint[3] = boundary[0]->boundary[start]->child[0]->vertex[1]; 00205 edge_midpoint[4] = boundary[0]->boundary[ii[start+step]]->child[0]->vertex[1]; 00206 edge_midpoint[5] = boundary[0]->boundary[ii[start+2*step]]->child[0]->vertex[1]; 00207 00208 surface_child[0][0] = boundary[0]->child[start]; 00209 surface_child[0][1] = boundary[0]->child[ii[start+step]]; 00210 surface_child[0][2] = boundary[0]->child[ii[start+2*step]]; 00211 surface_child[0][3] = boundary[0]->child[3]; 00212 00213 // for the second surface of the tetrahedron, we should judge its pose at first, too. 00214 // boundary[1]: the second surface 00215 // boundary[1]->vertex[i]: the i-th vertex of the triangle 00216 for (start = 0;start < 3; start ++) 00217 if (boundary[1]->vertex[start] == vertex[0]) 00218 break; 00219 Assert (start < 3,ExcInternalError()); 00220 if (boundary[1]->vertex[ii[start+1]] == vertex[2]) 00221 step = 1; 00222 else { 00223 Assert(boundary[1]->vertex[ii[start+1]] == vertex[3], ExcInternalError()); 00224 Assert(boundary[1]->vertex[ii[start+2]] == vertex[2], ExcInternalError()); 00225 step = 2; 00226 } 00227 Assert(edge_midpoint[3] == boundary[1]->boundary[start]->child[0]->vertex[1], ExcInternalError()); 00228 edge_midpoint[2] = boundary[1]->boundary[ii[start+step]]->child[0]->vertex[1]; 00229 edge_midpoint[1] = boundary[1]->boundary[ii[start+2*step]]->child[0]->vertex[1]; 00230 00231 surface_child[1][0] = boundary[1]->child[start]; 00232 surface_child[1][1] = boundary[1]->child[ii[start+step]]; 00233 surface_child[1][2] = boundary[1]->child[ii[start+2*step]]; 00234 surface_child[1][3] = boundary[1]->child[3]; 00235 00236 // for the third surface of the tetrahedron, we should judge its pose at first, too. 00237 // boundary[2]: the third surface 00238 // boundary[2]->vertex[i]: the i-th vertex of the triangle 00239 for (start = 0;start < 3; start ++) 00240 if (boundary[2]->vertex[start] == vertex[0]) 00241 break; 00242 Assert (start < 3,ExcInternalError()); 00243 if (boundary[2]->vertex[ii[start+1]] == vertex[1]) 00244 step = 1; 00245 else { 00246 Assert(boundary[2]->vertex[ii[start+1]] == vertex[3], ExcInternalError()); 00247 Assert(boundary[2]->vertex[ii[start+2]] == vertex[1], ExcInternalError()); 00248 step = 2; 00249 } 00250 Assert(edge_midpoint[4] == boundary[2]->boundary[start]->child[0]->vertex[1], ExcInternalError()); 00251 Assert(edge_midpoint[2] == boundary[2]->boundary[ii[start+step]]->child[0]->vertex[1], ExcInternalError()); 00252 edge_midpoint[0] = boundary[2]->boundary[ii[start+2*step]]->child[0]->vertex[1]; 00253 00254 surface_child[2][0] = boundary[2]->child[start]; 00255 surface_child[2][1] = boundary[2]->child[ii[start+step]]; 00256 surface_child[2][2] = boundary[2]->child[ii[start+2*step]]; 00257 surface_child[2][3] = boundary[2]->child[3]; 00258 00259 // for the fourth surface of the tetrahedron, we should judge its pose at first, too. 00260 // boundary[3]: the fourth surface 00261 // boundary[3]->vertex[i]: the i-th vertex of the triangle 00262 for (start = 0;start < 3; start ++) 00263 if (boundary[3]->vertex[start] == vertex[0]) 00264 break; 00265 Assert (start < 3,ExcInternalError()); 00266 if (boundary[3]->vertex[ii[start+1]] == vertex[1]) 00267 step = 1; 00268 else { 00269 Assert(boundary[3]->vertex[ii[start+1]] == vertex[2], ExcInternalError()); 00270 Assert(boundary[3]->vertex[ii[start+2]] == vertex[1], ExcInternalError()); 00271 step = 2; 00272 } 00273 Assert(edge_midpoint[5] == boundary[3]->boundary[start]->child[0]->vertex[1], ExcInternalError()); 00274 Assert(edge_midpoint[1] == boundary[3]->boundary[ii[start+step]]->child[0]->vertex[1], ExcInternalError()); 00275 Assert(edge_midpoint[0] == boundary[3]->boundary[ii[start+2*step]]->child[0]->vertex[1], ExcInternalError()); 00276 00277 surface_child[3][0] = boundary[3]->child[start]; 00278 surface_child[3][1] = boundary[3]->child[ii[start+step]]; 00279 surface_child[3][2] = boundary[3]->child[ii[start+2*step]]; 00280 surface_child[3][3] = boundary[3]->child[3]; 00281 } 00282 00283 // we add those geometries should be added for all refine model at first, which 00284 // include four triangles and four tetrahedrons at the four vertices. 00285 HGeometry<2,DOW> * triangle[8]; 00286 for (i = 0;i < 8;i ++) { 00287 triangle[i] = new HGeometry<2,DOW>(); 00288 Assert(triangle[i] != NULL, ExcOutOfMemory()); 00289 } 00290 triangle[0]->vertex[0] = edge_midpoint[0]; 00291 triangle[0]->vertex[1] = edge_midpoint[1]; 00292 triangle[0]->vertex[2] = edge_midpoint[2]; 00293 triangle[0]->boundary[0] = surface_child[1][0]->boundary[0]; 00294 triangle[0]->boundary[1] = surface_child[2][0]->boundary[0]; 00295 triangle[0]->boundary[2] = surface_child[3][0]->boundary[0]; 00296 triangle[0]->bmark = bmark; 00297 00298 triangle[1]->vertex[0] = edge_midpoint[4]; 00299 triangle[1]->vertex[1] = edge_midpoint[5]; 00300 triangle[1]->vertex[2] = edge_midpoint[0]; 00301 triangle[1]->boundary[0] = surface_child[3][1]->boundary[0]; 00302 triangle[1]->boundary[1] = surface_child[2][1]->boundary[0]; 00303 triangle[1]->boundary[2] = surface_child[0][0]->boundary[0]; 00304 triangle[1]->bmark = bmark; 00305 00306 triangle[2]->vertex[0] = edge_midpoint[5]; 00307 triangle[2]->vertex[1] = edge_midpoint[3]; 00308 triangle[2]->vertex[2] = edge_midpoint[1]; 00309 triangle[2]->boundary[0] = surface_child[1][1]->boundary[0]; 00310 triangle[2]->boundary[1] = surface_child[3][2]->boundary[0]; 00311 triangle[2]->boundary[2] = surface_child[0][1]->boundary[0]; 00312 triangle[2]->bmark = bmark; 00313 00314 triangle[3]->vertex[0] = edge_midpoint[2]; 00315 triangle[3]->vertex[1] = edge_midpoint[3]; 00316 triangle[3]->vertex[2] = edge_midpoint[4]; 00317 triangle[3]->boundary[0] = surface_child[0][2]->boundary[0]; 00318 triangle[3]->boundary[1] = surface_child[2][2]->boundary[0]; 00319 triangle[3]->boundary[2] = surface_child[1][2]->boundary[0]; 00320 triangle[3]->bmark = bmark; 00321 00322 // the four child tetrahedrons on the four verteics 00323 for (i = 0;i < 8;i ++) { 00324 child[i] = new HGeometry<3,DOW>(); 00325 Assert (child[i] != NULL, ExcOutOfMemory()); 00326 } 00327 00328 child[0]->parent = this; 00329 child[0]->vertex[0] = vertex[0]; 00330 child[0]->vertex[1] = edge_midpoint[0]; 00331 child[0]->vertex[2] = edge_midpoint[1]; 00332 child[0]->vertex[3] = edge_midpoint[2]; 00333 child[0]->boundary[0] = triangle[0]; 00334 child[0]->boundary[1] = surface_child[1][0]; 00335 child[0]->boundary[2] = surface_child[2][0]; 00336 child[0]->boundary[3] = surface_child[3][0]; 00337 child[0]->bmark = bmark; 00338 00339 child[1]->parent = this; 00340 child[1]->vertex[0] = vertex[1]; 00341 child[1]->vertex[1] = edge_midpoint[5]; 00342 child[1]->vertex[2] = edge_midpoint[0]; 00343 child[1]->vertex[3] = edge_midpoint[4]; 00344 child[1]->boundary[0] = triangle[1]; 00345 child[1]->boundary[1] = surface_child[2][1]; 00346 child[1]->boundary[2] = surface_child[0][0]; 00347 child[1]->boundary[3] = surface_child[3][1]; 00348 child[1]->bmark = bmark; 00349 00350 child[2]->parent = this; 00351 child[2]->vertex[0] = vertex[2]; 00352 child[2]->vertex[1] = edge_midpoint[1]; 00353 child[2]->vertex[2] = edge_midpoint[5]; 00354 child[2]->vertex[3] = edge_midpoint[3]; 00355 child[2]->boundary[0] = triangle[2]; 00356 child[2]->boundary[1] = surface_child[0][1]; 00357 child[2]->boundary[2] = surface_child[1][1]; 00358 child[2]->boundary[3] = surface_child[3][2]; 00359 child[2]->bmark = bmark; 00360 00361 child[3]->parent = this; 00362 child[3]->vertex[0] = vertex[3]; 00363 child[3]->vertex[1] = edge_midpoint[2]; 00364 child[3]->vertex[2] = edge_midpoint[3]; 00365 child[3]->vertex[3] = edge_midpoint[4]; 00366 child[3]->boundary[0] = triangle[3]; 00367 child[3]->boundary[1] = surface_child[0][2]; 00368 child[3]->boundary[2] = surface_child[2][2]; 00369 child[3]->boundary[3] = surface_child[1][2]; 00370 child[3]->bmark = bmark; 00371 00372 // choose the best refinement model. There are total 3 models to refine the 00373 // central part of the tetrahedron, between which the key difference is which 00374 // line is chosed as the main axis of geometry. To make the obtained tetrahedrons 00375 // more regular, we will choose the shortest diagonal line as the main axis. 00376 double l03, l14, l25; 00377 l03 = (*edge_midpoint[0] - *edge_midpoint[3]).length(); 00378 l14 = (*edge_midpoint[1] - *edge_midpoint[4]).length(); 00379 l25 = (*edge_midpoint[2] - *edge_midpoint[5]).length(); 00380 if (l03 <= l14) 00381 if (l03 <= l25) 00382 refine_model = REFINE_MODEL_03; 00383 else 00384 refine_model = REFINE_MODEL_25; 00385 else 00386 if (l14 <= l25) 00387 refine_model = REFINE_MODEL_14; 00388 else 00389 refine_model = REFINE_MODEL_25; 00390 00391 HGeometry<1,DOW> * axis; 00392 switch (refine_model) { 00393 00394 case REFINE_MODEL_03: 00395 // add the main axis at first 00396 axis = new HGeometry<1,DOW>(); 00397 Assert (axis != NULL, ExcOutOfMemory()); 00398 axis->vertex[0] = edge_midpoint[0]; 00399 axis->vertex[1] = edge_midpoint[3]; 00400 axis->bmark = bmark; 00401 00402 // add the four triangles 00403 triangle[4]->vertex[0] = edge_midpoint[4]; 00404 triangle[4]->vertex[1] = edge_midpoint[3]; 00405 triangle[4]->vertex[2] = edge_midpoint[0]; 00406 triangle[4]->boundary[0] = axis; 00407 triangle[4]->boundary[1] = surface_child[2][1]->boundary[0]; 00408 triangle[4]->boundary[2] = surface_child[0][2]->boundary[0]; 00409 triangle[4]->bmark = bmark; 00410 00411 triangle[5]->vertex[0] = edge_midpoint[2]; 00412 triangle[5]->vertex[1] = edge_midpoint[3]; 00413 triangle[5]->vertex[2] = edge_midpoint[0]; 00414 triangle[5]->boundary[0] = axis; 00415 triangle[5]->boundary[1] = surface_child[2][0]->boundary[0]; 00416 triangle[5]->boundary[2] = surface_child[1][2]->boundary[0]; 00417 triangle[5]->bmark = bmark; 00418 00419 triangle[6]->vertex[0] = edge_midpoint[1]; 00420 triangle[6]->vertex[1] = edge_midpoint[3]; 00421 triangle[6]->vertex[2] = edge_midpoint[0]; 00422 triangle[6]->boundary[0] = axis; 00423 triangle[6]->boundary[1] = surface_child[3][0]->boundary[0]; 00424 triangle[6]->boundary[2] = surface_child[1][1]->boundary[0]; 00425 triangle[6]->bmark = bmark; 00426 00427 triangle[7]->vertex[0] = edge_midpoint[5]; 00428 triangle[7]->vertex[1] = edge_midpoint[3]; 00429 triangle[7]->vertex[2] = edge_midpoint[0]; 00430 triangle[7]->boundary[0] = axis; 00431 triangle[7]->boundary[1] = surface_child[3][1]->boundary[0]; 00432 triangle[7]->boundary[2] = surface_child[0][1]->boundary[0]; 00433 triangle[7]->bmark = bmark; 00434 00435 // add the four tetrahedrons 00436 child[4]->parent = this; 00437 child[4]->vertex[0] = edge_midpoint[0]; 00438 child[4]->vertex[1] = edge_midpoint[3]; 00439 child[4]->vertex[2] = edge_midpoint[4]; 00440 child[4]->vertex[3] = edge_midpoint[5]; 00441 child[4]->boundary[0] = surface_child[0][3]; 00442 child[4]->boundary[1] = triangle[1]; 00443 child[4]->boundary[2] = triangle[7]; 00444 child[4]->boundary[3] = triangle[4]; 00445 child[4]->bmark = bmark; 00446 00447 child[5]->parent = this; 00448 child[5]->vertex[0] = edge_midpoint[0]; 00449 child[5]->vertex[1] = edge_midpoint[3]; 00450 child[5]->vertex[2] = edge_midpoint[1]; 00451 child[5]->vertex[3] = edge_midpoint[2]; 00452 child[5]->boundary[0] = surface_child[1][3]; 00453 child[5]->boundary[1] = triangle[0]; 00454 child[5]->boundary[2] = triangle[5]; 00455 child[5]->boundary[3] = triangle[6]; 00456 child[5]->bmark = bmark; 00457 00458 child[6]->parent = this; 00459 child[6]->vertex[0] = edge_midpoint[3]; 00460 child[6]->vertex[1] = edge_midpoint[0]; 00461 child[6]->vertex[2] = edge_midpoint[4]; 00462 child[6]->vertex[3] = edge_midpoint[2]; 00463 child[6]->boundary[0] = surface_child[2][3]; 00464 child[6]->boundary[1] = triangle[3]; 00465 child[6]->boundary[2] = triangle[5]; 00466 child[6]->boundary[3] = triangle[4]; 00467 child[6]->bmark = bmark; 00468 00469 child[7]->parent = this; 00470 child[7]->vertex[0] = edge_midpoint[3]; 00471 child[7]->vertex[1] = edge_midpoint[0]; 00472 child[7]->vertex[2] = edge_midpoint[1]; 00473 child[7]->vertex[3] = edge_midpoint[5]; 00474 child[7]->boundary[0] = surface_child[3][3]; 00475 child[7]->boundary[1] = triangle[2]; 00476 child[7]->boundary[2] = triangle[7]; 00477 child[7]->boundary[3] = triangle[6]; 00478 child[7]->bmark = bmark; 00479 00480 break; 00481 00482 case REFINE_MODEL_14: 00483 // add the main axis at first 00484 axis = new HGeometry<1,DOW>(); 00485 Assert (axis != NULL, ExcOutOfMemory()); 00486 axis->vertex[0] = edge_midpoint[1]; 00487 axis->vertex[1] = edge_midpoint[4]; 00488 axis->bmark = bmark; 00489 00490 // add the four triangles 00491 triangle[4]->vertex[0] = edge_midpoint[5]; 00492 triangle[4]->vertex[1] = edge_midpoint[4]; 00493 triangle[4]->vertex[2] = edge_midpoint[1]; 00494 triangle[4]->boundary[0] = axis; 00495 triangle[4]->boundary[1] = surface_child[3][2]->boundary[0]; 00496 triangle[4]->boundary[2] = surface_child[0][0]->boundary[0]; 00497 triangle[4]->bmark = bmark; 00498 00499 triangle[5]->vertex[0] = edge_midpoint[0]; 00500 triangle[5]->vertex[1] = edge_midpoint[4]; 00501 triangle[5]->vertex[2] = edge_midpoint[1]; 00502 triangle[5]->boundary[0] = axis; 00503 triangle[5]->boundary[1] = surface_child[3][0]->boundary[0]; 00504 triangle[5]->boundary[2] = surface_child[2][1]->boundary[0]; 00505 triangle[5]->bmark = bmark; 00506 00507 triangle[6]->vertex[0] = edge_midpoint[2]; 00508 triangle[6]->vertex[1] = edge_midpoint[4]; 00509 triangle[6]->vertex[2] = edge_midpoint[1]; 00510 triangle[6]->boundary[0] = axis; 00511 triangle[6]->boundary[1] = surface_child[1][0]->boundary[0]; 00512 triangle[6]->boundary[2] = surface_child[2][2]->boundary[0]; 00513 triangle[6]->bmark = bmark; 00514 00515 triangle[7]->vertex[0] = edge_midpoint[3]; 00516 triangle[7]->vertex[1] = edge_midpoint[4]; 00517 triangle[7]->vertex[2] = edge_midpoint[1]; 00518 triangle[7]->boundary[0] = axis; 00519 triangle[7]->boundary[1] = surface_child[1][1]->boundary[0]; 00520 triangle[7]->boundary[2] = surface_child[0][2]->boundary[0]; 00521 triangle[7]->bmark = bmark; 00522 00523 // add the four tetrahedrons 00524 child[4]->parent = this; 00525 child[4]->vertex[0] = edge_midpoint[1]; 00526 child[4]->vertex[1] = edge_midpoint[4]; 00527 child[4]->vertex[2] = edge_midpoint[5]; 00528 child[4]->vertex[3] = edge_midpoint[3]; 00529 child[4]->boundary[0] = surface_child[0][3]; 00530 child[4]->boundary[1] = triangle[2]; 00531 child[4]->boundary[2] = triangle[7]; 00532 child[4]->boundary[3] = triangle[4]; 00533 child[4]->bmark = bmark; 00534 00535 child[5]->parent = this; 00536 child[5]->vertex[0] = edge_midpoint[4]; 00537 child[5]->vertex[1] = edge_midpoint[1]; 00538 child[5]->vertex[2] = edge_midpoint[2]; 00539 child[5]->vertex[3] = edge_midpoint[3]; 00540 child[5]->boundary[0] = surface_child[1][3]; 00541 child[5]->boundary[1] = triangle[3]; 00542 child[5]->boundary[2] = triangle[7]; 00543 child[5]->boundary[3] = triangle[6]; 00544 child[5]->bmark = bmark; 00545 00546 child[6]->parent = this; 00547 child[6]->vertex[0] = edge_midpoint[1]; 00548 child[6]->vertex[1] = edge_midpoint[4]; 00549 child[6]->vertex[2] = edge_midpoint[2]; 00550 child[6]->vertex[3] = edge_midpoint[0]; 00551 child[6]->boundary[0] = surface_child[2][3]; 00552 child[6]->boundary[1] = triangle[0]; 00553 child[6]->boundary[2] = triangle[5]; 00554 child[6]->boundary[3] = triangle[6]; 00555 child[6]->bmark = bmark; 00556 00557 child[7]->parent = this; 00558 child[7]->vertex[0] = edge_midpoint[4]; 00559 child[7]->vertex[1] = edge_midpoint[1]; 00560 child[7]->vertex[2] = edge_midpoint[5]; 00561 child[7]->vertex[3] = edge_midpoint[0]; 00562 child[7]->boundary[0] = surface_child[3][3]; 00563 child[7]->boundary[1] = triangle[1]; 00564 child[7]->boundary[2] = triangle[5]; 00565 child[7]->boundary[3] = triangle[4]; 00566 child[7]->bmark = bmark; 00567 00568 break; 00569 00570 case REFINE_MODEL_25: 00571 // add the main axis at first 00572 axis = new HGeometry<1,DOW>(); 00573 Assert (axis != NULL, ExcOutOfMemory()); 00574 axis->vertex[0] = edge_midpoint[2]; 00575 axis->vertex[1] = edge_midpoint[5]; 00576 axis->bmark = bmark; 00577 00578 // add the four triangles 00579 triangle[4]->vertex[0] = edge_midpoint[4]; 00580 triangle[4]->vertex[1] = edge_midpoint[5]; 00581 triangle[4]->vertex[2] = edge_midpoint[2]; 00582 triangle[4]->boundary[0] = axis; 00583 triangle[4]->boundary[1] = surface_child[2][2]->boundary[0]; 00584 triangle[4]->boundary[2] = surface_child[0][0]->boundary[0]; 00585 triangle[4]->bmark = bmark; 00586 00587 triangle[5]->vertex[0] = edge_midpoint[0]; 00588 triangle[5]->vertex[1] = edge_midpoint[5]; 00589 triangle[5]->vertex[2] = edge_midpoint[2]; 00590 triangle[5]->boundary[0] = axis; 00591 triangle[5]->boundary[1] = surface_child[2][0]->boundary[0]; 00592 triangle[5]->boundary[2] = surface_child[3][1]->boundary[0]; 00593 triangle[5]->bmark = bmark; 00594 00595 triangle[6]->vertex[0] = edge_midpoint[1]; 00596 triangle[6]->vertex[1] = edge_midpoint[5]; 00597 triangle[6]->vertex[2] = edge_midpoint[2]; 00598 triangle[6]->boundary[0] = axis; 00599 triangle[6]->boundary[1] = surface_child[1][0]->boundary[0]; 00600 triangle[6]->boundary[2] = surface_child[3][2]->boundary[0]; 00601 triangle[6]->bmark = bmark; 00602 00603 triangle[7]->vertex[0] = edge_midpoint[3]; 00604 triangle[7]->vertex[1] = edge_midpoint[5]; 00605 triangle[7]->vertex[2] = edge_midpoint[2]; 00606 triangle[7]->boundary[0] = axis; 00607 triangle[7]->boundary[1] = surface_child[1][2]->boundary[0]; 00608 triangle[7]->boundary[2] = surface_child[0][1]->boundary[0]; 00609 triangle[7]->bmark = bmark; 00610 00611 // add the four tetrahedrons 00612 child[4]->parent = this; 00613 child[4]->vertex[0] = edge_midpoint[2]; 00614 child[4]->vertex[1] = edge_midpoint[5]; 00615 child[4]->vertex[2] = edge_midpoint[3]; 00616 child[4]->vertex[3] = edge_midpoint[4]; 00617 child[4]->boundary[0] = surface_child[0][3]; 00618 child[4]->boundary[1] = triangle[3]; 00619 child[4]->boundary[2] = triangle[4]; 00620 child[4]->boundary[3] = triangle[7]; 00621 child[4]->bmark = bmark; 00622 00623 child[5]->parent = this; 00624 child[5]->vertex[0] = edge_midpoint[5]; 00625 child[5]->vertex[1] = edge_midpoint[2]; 00626 child[5]->vertex[2] = edge_midpoint[3]; 00627 child[5]->vertex[3] = edge_midpoint[1]; 00628 child[5]->boundary[0] = surface_child[1][3]; 00629 child[5]->boundary[1] = triangle[2]; 00630 child[5]->boundary[2] = triangle[6]; 00631 child[5]->boundary[3] = triangle[7]; 00632 child[5]->bmark = bmark; 00633 00634 child[6]->parent = this; 00635 child[6]->vertex[0] = edge_midpoint[5]; 00636 child[6]->vertex[1] = edge_midpoint[2]; 00637 child[6]->vertex[2] = edge_midpoint[0]; 00638 child[6]->vertex[3] = edge_midpoint[4]; 00639 child[6]->boundary[0] = surface_child[2][3]; 00640 child[6]->boundary[1] = triangle[1]; 00641 child[6]->boundary[2] = triangle[4]; 00642 child[6]->boundary[3] = triangle[5]; 00643 child[6]->bmark = bmark; 00644 00645 child[7]->parent = this; 00646 child[7]->vertex[0] = edge_midpoint[2]; 00647 child[7]->vertex[1] = edge_midpoint[5]; 00648 child[7]->vertex[2] = edge_midpoint[0]; 00649 child[7]->vertex[3] = edge_midpoint[1]; 00650 child[7]->boundary[0] = surface_child[3][3]; 00651 child[7]->boundary[1] = triangle[0]; 00652 child[7]->boundary[2] = triangle[6]; 00653 child[7]->boundary[3] = triangle[5]; 00654 child[7]->bmark = bmark; 00655 00656 break; 00657 00658 default: 00659 Assert(false, ExcInternalError()); 00660 } 00661 } 00662 00663 00664 template <> 00665 void HGeometry<1,DOW>::checkIntegrity() const 00666 { 00667 if (!isRefined()) return; 00668 for (int i = 0;i < n_child;i ++) { 00669 child[i]->checkIntegrity(); 00670 } 00671 } 00672 00673 00674 00675 template <> 00676 void HGeometry<2,DOW>::checkIntegrity() const 00677 { 00678 int i, j, k; 00679 for (i = 0;i < n_boundary;i ++) { 00680 HGeometry<1,DOW> * b = boundary[i]; 00681 for (j = 0;j < b->n_vertex;j ++) { 00682 for (k = 0;k < n_vertex;k ++) { 00683 if (k == i) continue; 00684 if (b->vertex[j] == vertex[k]) 00685 break; 00686 } 00687 Assert(k < n_vertex, ExcInternalError()); 00688 } 00689 } 00690 if (!isRefined()) return; 00691 for (int i = 0;i < n_child;i ++) { 00692 child[i]->checkIntegrity(); 00693 } 00694 } 00695 00696 00697 00698 template <> 00699 void HGeometry<3,DOW>::checkIntegrity() const 00700 { 00701 int i, j, k; 00702 for (i = 0;i < n_boundary;i ++) { 00703 HGeometry<2,DOW> * b = boundary[i]; 00704 for (j = 0;j < b->n_vertex;j ++) { 00705 for (k = 0;k < n_vertex;k ++) { 00706 if (k == i) continue; 00707 if (b->vertex[j] == vertex[k]) 00708 break; 00709 } 00710 Assert(k < n_vertex, ExcInternalError()); 00711 } 00712 } 00713 if (!isRefined()) return; 00714 for (int i = 0;i < n_child;i ++) { 00715 child[i]->checkIntegrity(); 00716 } 00717 } 00718 00719 00720 00721 template <> 00722 bool HGeometry<1,DOW>::isIncludePoint(const afepack::Point<DOW>& p) const 00723 { 00724 std::cerr << "warning: HGeometry<1,DOW>::isIncludePoint not implemented, return true" << std::endl; 00725 return true; 00726 } 00727 00728 00729 template <> 00730 bool HGeometry<2,DOW>::isIncludePoint(const afepack::Point<DOW>& p) const 00731 { 00732 double lambda[3]; 00733 const std::vector<vertex_t *>& v = vertex; 00734 double volume = (((*v[1])[0] - (*v[0])[0])*((*v[2])[1] - (*v[0])[1]) - 00735 ((*v[1])[1] - (*v[0])[1])*((*v[2])[0] - (*v[0])[0])); 00736 double mzero = -1.0e-08; 00737 lambda[0] = (((*v[1])[0] - p[0])*((*v[2])[1] - p[1]) - 00738 ((*v[1])[1] - p[1])*((*v[2])[0] - p[0]))/volume; 00739 lambda[1] = (((*v[2])[0] - p[0])*((*v[0])[1] - p[1]) - 00740 ((*v[2])[1] - p[1])*((*v[0])[0] - p[0]))/volume; 00741 lambda[2] = (((*v[0])[0] - p[0])*((*v[1])[1] - p[1]) - 00742 ((*v[0])[1] - p[1])*((*v[1])[0] - p[0]))/volume; 00743 return (lambda[0] >= mzero && lambda[1] >= mzero && lambda[2] >= mzero); 00744 } 00745 00746 00747 00748 template <> 00749 bool HGeometry<3,DOW>::isIncludePoint(const afepack::Point<DOW>& p) const 00750 { 00751 double ZERO = 1.0e-06; 00752 const std::vector<vertex_t *>& v = vertex; 00753 double A[4] = { 00754 tetrahedron_volume( p, *v[1], *v[2], *v[3]), 00755 tetrahedron_volume(*v[0], p, *v[2], *v[3]), 00756 tetrahedron_volume(*v[0], *v[1], p, *v[3]), 00757 tetrahedron_volume(*v[0], *v[1], *v[2], p) 00758 }; 00759 double AZERO = ZERO*(A[0] + A[1] + A[2] + A[3]); 00760 return ((A[0] >= -AZERO) && (A[1] >= -AZERO) && 00761 (A[2] >= -AZERO) && (A[3] >= -AZERO)); 00762 } 00763 00764 00765 00766 template <> 00767 void HElement<2, DOW>::refine() 00768 { 00769 if (isRefined()) return; 00770 00771 h_element->refine(); 00772 00773 // construct its own child 00774 for (int i = 0;i < n_child;i ++) { 00775 child[i] = new HElement<2, DOW>(); 00776 Assert (child[i] != NULL, ExcOutOfMemory()); 00777 child[i]->h_element = h_element->child[i]; 00778 child[i]->parent = this; 00779 } 00780 } 00781 00782 00783 00784 template <> 00785 void HElement<3, DOW>::refine() 00786 { 00787 if (isRefined()) return; 00788 00789 // refine the h geometry tree at first 00790 h_element->refine(); 00791 00792 // construct its own child 00793 for (int i = 0;i < n_child;i ++) { 00794 child[i] = new HElement<3, DOW>(); 00795 Assert (child[i] != NULL, ExcOutOfMemory()); 00796 child[i]->h_element = h_element->child[i]; 00797 child[i]->parent = this; 00798 } 00799 } 00800 00801 00802 00803 template <> 00804 std::ostream& operator<<(std::ostream& os, const HGeometry<1,DOW>& geometry) 00805 { 00806 os << (*geometry.vertex[0])[0] << "\t" 00807 << (*geometry.vertex[1])[0] << "\t" 00808 << (*geometry.vertex[0])[1] << "\t" 00809 << (*geometry.vertex[1])[1] << "\n"; 00810 return os; 00811 } 00812 00813 00814 00815 template <> 00816 std::ostream& operator<<(std::ostream& os, const HGeometry<0,DOW>& geometry) 00817 { 00818 // for (int i = 0;i < DOW;i ++) 00819 // os << geometry[i] << "\t"; 00820 return os; 00821 } 00822 00823 00824 00825 template <> 00826 void IrregularMesh<2, DOW>::regularize(bool renumerate) 00827 { 00828 if (regular_mesh != NULL) 00829 delete regular_mesh; 00830 std::cerr << "Generating regular mesh from the semiregular mesh ..." << std::flush; 00831 int i, j, k, n_node, n_side, n_element; 00832 int ii[] = {0, 1, 2, 0, 1, 2, 0, 1, 2}; 00833 00834 std::cerr << "\n\tpreparing data ..." << std::flush; 00835 Tools tools; 00836 ActiveElementIterator<2, DOW> 00837 the_ele = beginActiveElement(), 00838 end_ele = endActiveElement(); 00839 // ΪԼ 00840 n_element = 0; 00841 for (;the_ele != end_ele;++ the_ele) { 00842 HGeometry<2, DOW> *& h_element = the_ele->h_element; 00847 for (i = 0;i < h_element->n_vertex;i ++) { 00848 tools.setGeometryActive(*(h_element->vertex[i])); 00849 } 00855 int n_refined_edge = 0; 00856 for (i = 0;i < the_ele->n_boundary;i ++) { 00857 HGeometry<1, DOW>& side = *(h_element->boundary[i]); 00862 if (tools.isRefined(side)) { 00863 tools.setGeometryActive(*side.child[0]->vertex[1]); 00864 tools.setGeometryActive(*side.child[0]); 00865 tools.setGeometryActive(*side.child[1]); 00866 tools.setGeometryInactive(side); 00867 n_refined_edge += 1; 00868 } else { 00869 tools.setGeometryActive(side); 00870 } 00871 } 00872 assert (n_refined_edge <= 1); 00873 assert (tools.isGeometryUsed(*h_element)); 00874 tools.setParentInactive(*h_element); 00875 n_element += 1; 00876 } 00877 00878 // м 00879 n_node = 0; 00880 n_side = 0; 00881 n_element = 0; 00882 the_ele = beginActiveElement(); 00883 for (;the_ele != end_ele;++ the_ele) { 00884 HGeometry<2, DOW> *& h_element = the_ele->h_element; 00889 for (i = 0;i < the_ele->n_vertex;i ++) { 00890 HGeometry<0, DOW>& node = *(h_element->vertex[i]); 00891 if (tools.isGeometryActive(node)) { 00892 node.index = n_node ++; 00893 } 00894 } 00899 int n_refined_edge = 0; 00900 for (i = 0;i < the_ele->n_boundary;i ++) { 00901 HGeometry<1, DOW>& side = *(h_element->boundary[i]); 00902 if (tools.isGeometryActive(side)) { 00903 side.index = n_side ++; 00904 } else if (tools.isGeometryInactive(side)) { 00905 HGeometry<0, DOW>& mp = *(side.child[0]->vertex[1]); 00906 if (tools.isGeometryActive(mp)) { 00907 mp.index = n_node ++; 00908 } 00909 for (j = 0;j < side.n_child;++ j) { 00910 HGeometry<1, DOW> *& chd = side.child[j]; 00911 if (tools.isGeometryActive(*chd)) { 00912 chd->index = n_side ++; 00913 } 00914 } 00915 n_refined_edge += 1; 00916 } 00917 } 00918 assert (n_refined_edge <= 1); 00919 assert (tools.isGeometryActive(*h_element)); 00920 h_element->index = n_element ++; 00921 the_ele->index = h_element->index; 00922 } 00923 00924 std::cerr << "\n\tbuilding the regular mesh ..." << std::flush; 00925 // 00926 regular_mesh = new RegularMesh<2, DOW>(this); 00927 Assert (regular_mesh != NULL, ExcOutOfMemory()); 00928 regular_mesh->point().resize(n_node); 00929 regular_mesh->geometry(0).resize(n_node); 00930 regular_mesh->geometry(1).resize(n_side); 00931 regular_mesh->geometry(2).resize(n_element); 00932 00933 #ifdef __SERIALIZATION__ 00934 std::vector<std::vector<HGeometryBase*> >& h_geometry = regular_mesh->h_geometry(); 00935 h_geometry.clear(); 00936 h_geometry.resize(3); 00937 h_geometry[0].resize(n_node); 00938 h_geometry[1].resize(n_side); 00939 h_geometry[2].resize(n_element); 00940 #endif 00941 00942 int n_triangle = 0, n_twin_triangle = 0; 00943 the_ele = beginActiveElement(); 00944 for (;the_ele != end_ele;++ the_ele) { 00945 HGeometry<2, DOW> * h_element = the_ele->h_element; 00946 // ȼεĶ 00947 for (i = 0;i < the_ele->n_vertex;i ++) { 00948 HGeometry<0, DOW>& vtx = *(h_element->vertex[i]); 00949 j = vtx.index; 00950 GeometryBM& pnt = regular_mesh->geometry(0,j); 00951 tools.regularize_add_node(vtx, 00952 pnt, 00953 regular_mesh->point(j)); 00954 #ifdef __SERIALIZATION__ 00955 h_geometry[0][j] = &vtx; 00956 #endif 00957 } 00958 00959 // Ȼε 00960 int m_refined_edge = -1; 00961 for (i = 0;i < the_ele->n_boundary;i ++) { 00962 HGeometry<1, DOW>& bnd = *(h_element->boundary[i]); 00963 if (tools.isGeometryInactive(bnd)) { // ϸֵı 00964 assert (m_refined_edge == -1); 00965 m_refined_edge = i; 00966 00968 HGeometry<0, DOW>& mp = *(bnd.child[0]->vertex[1]); 00969 j = mp.index; 00970 tools.regularize_add_node(mp, 00971 regular_mesh->geometry(0,j), 00972 regular_mesh->point(j)); 00973 #ifdef __SERIALIZATION__ 00974 h_geometry[0][j] = ∓ 00975 #endif 00976 00977 HGeometry<1, DOW>& bnd0 = *(bnd.child[0]); 00978 j = bnd0.index; 00979 tools.regularize_add_side(bnd0, 00980 regular_mesh->geometry(1,j)); 00981 #ifdef __SERIALIZATION__ 00982 h_geometry[1][j] = &bnd0; 00983 #endif 00984 00985 HGeometry<1, DOW>& bnd1 = *(bnd.child[1]); 00986 j = bnd1.index; 00987 tools.regularize_add_side(bnd1, 00988 regular_mesh->geometry(1,j)); 00989 #ifdef __SERIALIZATION__ 00990 h_geometry[1][j] = &bnd1; 00991 #endif 00992 } else { 00993 00994 j = bnd.index; 00995 tools.regularize_add_side(bnd, 00996 regular_mesh->geometry(1,j)); 00997 #ifdef __SERIALIZATION__ 00998 h_geometry[1][j] = &bnd; 00999 #endif 01000 } 01001 } 01002 01003 // Ϊκ˫ 01004 j = h_element->index; 01005 assert (j == the_ele->index); 01006 GeometryBM& element = regular_mesh->geometry(2,j); 01007 Assert(element.index() == -1, ExcInternalError()); 01008 #ifdef __SERIALIZATION__ 01009 h_geometry[2][j] = h_element; 01010 #endif 01011 element.boundaryMark() = h_element->bmark; 01012 if (m_refined_edge == -1) { 01013 n_triangle ++; 01014 tools.regularize_add_triangle(*h_element, element); 01015 } else { // ˫ε 01016 n_twin_triangle ++; 01017 tools.regularize_add_twin_triangle(*h_element, element, m_refined_edge); 01018 } 01019 } 01020 01021 geometry_tree->unlock(); 01022 std::cerr << "\n\tnodes: " << n_node 01023 << "; sides: " << n_side 01024 << "; elements: " << n_element 01025 << " (" << n_triangle 01026 << ", " << n_twin_triangle 01027 << ")" << std::endl; 01028 01029 if (renumerate) renumerateElement(); 01030 } 01031 01032 01033 01034 template <> 01035 void IrregularMesh<3, DOW>::regularize(bool renumerate) 01036 { 01037 if (regular_mesh != NULL) 01038 delete regular_mesh; 01039 std::cerr << "Generating regular mesh from the semiregular mesh ..." << std::flush; 01040 int ii[] = {0, 1, 2, 0, 1, 2, 0, 1, 2}; 01041 Tools tools; 01042 01043 std::cerr << "\n\tpreparing data ..." << std::flush; 01044 ActiveElementIterator<3, DOW> 01045 the_ele = beginActiveElement(), 01046 end_ele = endActiveElement(); 01047 // Ϊм 01048 for (;the_ele != end_ele;++ the_ele) { 01049 HGeometry<3, DOW> * h_element = the_ele->h_element; 01050 for (int i = 0;i < the_ele->n_vertex;++ i) { 01051 tools.setGeometryActive(*(h_element->vertex[i])); 01052 } 01053 for (int i = 0;i < h_element->n_boundary;++ i) { 01054 HGeometry<2, DOW> * surface = h_element->boundary[i]; 01061 if (tools.isRefined(*surface)) { 01062 01063 for (int j = 0;j < surface->child[3]->n_boundary;++ j) { 01064 assert ((! tools.isRefined(*surface->child[3]->boundary[j]))); 01065 tools.setGeometryActive(*surface->child[3]->boundary[j]); 01066 } 01067 01069 for (int j = 0;j < surface->n_child;++ j) { 01070 assert ((! tools.isRefined(*surface->child[j]))); 01071 tools.setGeometryActive(*surface->child[j]); 01072 } 01073 01075 tools.setGeometryInactive(*surface); 01076 } else { 01077 tools.setGeometryActive(*surface); 01078 } 01079 int n_refined_edge = 0; 01080 for (int j = 0;j < surface->n_boundary;++ j) { 01081 HGeometry<1, DOW> * edge = surface->boundary[j]; 01082 if (tools.isRefined(*edge)) { 01083 tools.setGeometryActive(*edge->child[0]->vertex[1]); 01084 01085 for (int k = 0;k < edge->n_child;++ k) { 01086 assert ((! tools.isRefined(*edge->child[k]))); 01087 tools.setGeometryActive(*edge->child[k]); 01088 } 01089 01090 tools.setGeometryInactive(*edge); 01091 n_refined_edge += 1; 01092 } else { 01093 tools.setGeometryActive(*edge); 01094 } 01095 } 01096 if (!( (n_refined_edge == 0) || 01097 (n_refined_edge == 1) || 01098 ((n_refined_edge == 3) && 01099 tools.isRefined(*surface)))) { 01100 abort(); 01101 } 01102 } 01103 Assert(tools.isGeometryActive(*h_element), ExcInternalError()); 01104 if (! tools.isGeometryActive(*h_element)) abort(); 01105 tools.setParentInactive(*h_element); 01106 } 01107 01108 // м 01109 int n_node = 0; 01110 int n_edge = 0; 01111 int n_surface = 0; 01112 int n_element = 0; 01113 the_ele = beginActiveElement(); 01114 for (;the_ele != end_ele;++ the_ele) { 01115 HGeometry<3, DOW> * h_element = the_ele->h_element; 01116 // Ŷ 01117 for (int i = 0;i < the_ele->n_vertex;++ i) { 01118 HGeometry<0, DOW> * node = h_element->vertex[i]; 01119 if (tools.isGeometryActive(*node)) { 01120 node->index = n_node ++; 01121 } 01122 } 01123 01124 // 01125 for (int i = 0;i < the_ele->n_boundary;++ i) { 01126 HGeometry<2, DOW> * surface = h_element->boundary[i]; 01127 if (tools.isGeometryActive(*surface)) { 01128 surface->index = n_surface ++; 01129 01130 for (int j = 0;j < surface->n_boundary;++ j) { 01131 HGeometry<1, DOW> * edge = surface->boundary[j]; 01132 if (tools.isGeometryActive(*edge)) { 01133 edge->index = n_edge ++; 01134 } else if (tools.isGeometryInactive(*edge)) { 01143 for (int k = 0;k < edge->n_child;++ k) { 01144 HGeometry<1, DOW> * chd = edge->child[k]; 01145 if (tools.isGeometryActive(*chd)) { 01146 chd->index = n_edge ++; 01147 } 01148 } 01149 HGeometry<0,DOW> * vtx = edge->child[0]->vertex[1]; 01150 if (tools.isGeometryActive(*vtx)) { 01151 vtx->index = n_node ++; 01152 } 01153 } 01154 } 01155 } else if (tools.isGeometryInactive(*surface)) { 01156 for (int j = 0;j < surface->n_child;++ j) { 01157 HGeometry<2,DOW> * chd = surface->child[j]; 01158 if (tools.isGeometryActive(*chd)) { 01159 chd->index = n_surface ++; 01160 } 01161 for (int k = 0;k < chd->n_vertex;++ k) { 01162 HGeometry<0,DOW> * vtx = chd->vertex[k]; 01163 if (tools.isGeometryActive(*vtx)) { 01164 vtx->index = n_node ++; 01165 } 01166 } 01167 for (int k = 0;k < chd->n_boundary;++ k) { 01168 HGeometry<1,DOW> * bnd = chd->boundary[k]; 01169 if (tools.isGeometryActive(*bnd)) { 01170 bnd->index = n_edge ++; 01171 } 01172 } 01173 } 01174 } 01175 } 01176 Assert(tools.isGeometryActive(*h_element), ExcInternalError()); 01177 if (! tools.isGeometryActive(*h_element)) abort(); 01178 h_element->index = n_element ++; 01179 the_ele->index = h_element->index; 01180 } 01181 01182 std::cerr << "\n\tbuilding the regular mesh ..." << std::flush; 01183 // 01184 regular_mesh = new RegularMesh<3, DOW>(this); 01185 regular_mesh->point().resize(n_node); 01186 regular_mesh->geometry(0).resize(n_node); 01187 regular_mesh->geometry(1).resize(n_edge); 01188 regular_mesh->geometry(2).resize(n_surface); 01189 regular_mesh->geometry(3).resize(n_element); 01190 01191 #ifdef __SERIALIZATION__ 01192 std::vector<std::vector<HGeometryBase*> >& h_geometry = regular_mesh->h_geometry(); 01193 h_geometry.clear(); 01194 h_geometry.resize(4); 01195 h_geometry[0].resize(n_node); 01196 h_geometry[1].resize(n_edge); 01197 h_geometry[2].resize(n_surface); 01198 h_geometry[3].resize(n_element); 01199 #endif 01200 01201 int n_tetrahedron = 0; 01202 int n_twin_tetrahedron = 0; 01203 int n_four_tetrahedron = 0; 01204 the_ele = beginActiveElement(); 01205 for (;the_ele != end_ele;++ the_ele) { 01206 01207 HGeometry<3, DOW> * h_element = the_ele->h_element; 01208 // ȼ붥 01209 for (int i = 0;i < h_element->n_vertex;++ i) { 01210 HGeometry<0,DOW> * vtx = h_element->vertex[i]; 01211 int j = vtx->index; 01212 tools.regularize_add_node(*vtx, 01213 regular_mesh->geometry(0, j), 01214 regular_mesh->point(j)); 01215 #ifdef __SERIALIZATION__ 01216 h_geometry[0][j] = vtx; 01217 #endif 01218 } 01219 01220 // ȻԪı 01221 int n_twin_triangle_surface = 0, twin_triangle_surface[3]; 01222 int n_nonactive_neighbour = 0, nonactive_neighbour; 01223 HGeometry<1, DOW> * refined_edge = NULL; 01224 for (int i = 0;i < the_ele->n_boundary;++ i) { 01225 HGeometry<2, DOW> * bnd = h_element->boundary[i]; 01226 if (tools.isGeometryIndexed(*bnd)) { // һactive 01233 int k = 0; 01234 for (;k < bnd->n_boundary;k ++) { 01235 HGeometry<1,DOW> * edge = bnd->boundary[k]; 01236 if (tools.isGeometryInactive(*edge)) { 01237 refined_edge = edge; break; 01238 } 01239 } 01240 if (k < bnd->n_boundary) { 01241 twin_triangle_surface[n_twin_triangle_surface ++] = i; 01242 } 01243 01244 int j = bnd->index; 01245 GeometryBM& surface = regular_mesh->geometry(2,j); 01246 if (surface.index() != -1) continue; // 洦ˣͼһ 01247 01248 if (k == bnd->n_boundary) { // 01249 for (int l = 0;l < bnd->n_boundary;l ++) { // ȼ 01250 HGeometry<1,DOW> * edge = bnd->boundary[l]; 01251 int m = edge->index; 01252 tools.regularize_add_side(*edge, regular_mesh->geometry(1,m)); 01253 #ifdef __SERIALIZATION__ 01254 h_geometry[1][m] = edge; 01255 #endif 01256 } 01257 // ټ 01258 tools.regularize_add_triangle(*bnd, surface); 01259 #ifdef __SERIALIZATION__ 01260 h_geometry[2][j] = bnd; 01261 #endif 01262 } else { // ˫Σǵ k ߱ϸ 01266 for (int m = 1;m < 3;++ m) { 01267 HGeometry<1,DOW> * edge = bnd->boundary[ii[k+m]]; 01268 int n = edge->index; 01269 tools.regularize_add_side(*edge, regular_mesh->geometry(1,n)); 01270 #ifdef __SERIALIZATION__ 01271 h_geometry[1][n] = edge; 01272 #endif 01273 } 01274 01275 { 01276 for (int m = 0;m < refined_edge->n_child;++ m) { 01277 HGeometry<1,DOW> * edge = refined_edge->child[m]; 01278 int n = edge->index; 01279 tools.regularize_add_side(*edge, regular_mesh->geometry(1,n)); 01280 #ifdef __SERIALIZATION__ 01281 h_geometry[1][n] = edge; 01282 #endif 01283 } 01284 } 01285 01286 { 01287 HGeometry<0,DOW> * vtx = refined_edge->child[0]->vertex[1]; 01288 int n = vtx->index; 01289 tools.regularize_add_node(*vtx, 01290 regular_mesh->geometry(0, n), 01291 regular_mesh->point(n)); 01292 #ifdef __SERIALIZATION__ 01293 h_geometry[0][n] = vtx; 01294 #endif 01295 } 01296 01297 // Ȼ˫ 01298 tools.regularize_add_twin_triangle(*bnd, surface, k); 01299 #ifdef __SERIALIZATION__ 01300 h_geometry[2][j] = bnd; 01301 #endif 01302 } 01303 } else { 01307 n_nonactive_neighbour ++; 01308 Assert(n_nonactive_neighbour == 1, ExcInternalError()); 01309 nonactive_neighbour = i; // Ϳ 01310 01311 for (int j = 0;j < bnd->n_child;++ j) { 01312 HGeometry<2,DOW> * chd = bnd->child[j]; 01313 { 01314 int n = chd->index; 01315 tools.regularize_add_triangle(*chd, regular_mesh->geometry(2,n)); 01316 #ifdef __SERIALIZATION__ 01317 h_geometry[2][n] = chd; 01318 #endif 01319 } 01320 01321 for (int k = 0;k < chd->n_vertex;++ k) { 01322 HGeometry<0,DOW> * vtx = chd->vertex[k]; 01323 int n = vtx->index; 01324 tools.regularize_add_node(*vtx, 01325 regular_mesh->geometry(0, n), 01326 regular_mesh->point(n)); 01327 #ifdef __SERIALIZATION__ 01328 h_geometry[0][n] = vtx; 01329 #endif 01330 } 01331 01332 for (int k = 0;k < chd->n_boundary;++ k) { 01333 HGeometry<1,DOW> * edge = chd->boundary[k]; 01334 int n = edge->index; 01335 tools.regularize_add_side(*edge, 01336 regular_mesh->geometry(1,n)); 01337 #ifdef __SERIALIZATION__ 01338 h_geometry[1][n] = edge; 01339 #endif 01340 } 01341 } 01342 01343 } 01344 } 01345 01346 // 01347 int j = h_element->index; 01348 GeometryBM& tetra = regular_mesh->geometry(3,j); 01349 tetra.index() = j; 01350 #ifdef __SERIALIZATION__ 01351 h_geometry[3][j] = h_element; 01352 #endif 01353 if (n_nonactive_neighbour == 1) { // İ̥ 01354 Assert((n_twin_triangle_surface == 3), ExcInternalError()); 01355 01356 int tt[4][4] = {{0, 1, 2, 3}, {1, 2, 0, 3}, {2, 3, 0, 1}, {3, 0, 2, 1}}; 01357 int * t = &tt[nonactive_neighbour][0]; 01358 HGeometry<2,DOW> *& refined_face = h_element->boundary[t[0]]; 01359 int q[3]; 01360 for (int k = 0;k < 3;++ k) { 01361 for (int l = 0;l < 3;++ l) { 01362 if (h_element->vertex[t[k+1]] == refined_face->vertex[l]) { 01363 q[k] = l; break; 01364 } 01365 } 01366 } 01367 01368 tetra.vertex().resize(7); 01369 for (int k = 0;k < 4;++ k) { 01370 tetra.vertex(k) = h_element->vertex[t[k]]->index; 01371 } 01372 for (int k = 0;k < 3;++ k) { 01373 tetra.vertex(k+4) = refined_face->boundary[q[k]]->child[0]->vertex[1]->index; 01374 } 01375 01376 tetra.boundary().resize(7); 01377 tetra.boundary(0) = refined_face->child[3]->index; 01378 for (int k = 0;k < 3;++ k) { 01379 tetra.boundary(k+1) = refined_face->child[q[k]]->index; 01380 } 01381 for (int k = 0;k < 3;++ k) { 01382 tetra.boundary(k+4) = h_element->boundary[t[k+1]]->index; 01383 } 01384 01385 tetra.boundaryMark() = h_element->bmark; 01386 n_four_tetrahedron ++; 01387 } else if (n_twin_triangle_surface > 0) { // ˫̥ 01388 Assert (n_twin_triangle_surface == 2, ExcInternalError()); 01389 01390 // жϸ 01391 int rei = (1<<twin_triangle_surface[0]) + (1<<twin_triangle_surface[1]); 01392 int tt[16] = { 01393 -1, -1, -1, 3, -1, 4, 2, -1, 01394 -1, 5, 1, -1, 0, -1, -1, -1, 01395 }; 01396 int refined_edge_idx = tt[rei]; 01397 Assert ((refined_edge_idx >= 0), ExcInternalError()); 01398 01399 tetra.vertex().resize(5); 01400 tetra.boundary().resize(4); 01401 int kk[6][4] = { 01402 {3, 1, 0, 2}, {1, 2, 0, 3}, {1, 0, 3, 2}, 01403 {0, 2, 3, 1}, {0, 3, 1, 2}, {0, 1, 2, 3} 01404 }; 01405 for (int k = 0;k < 4;++ k) { 01406 int k0 = kk[refined_edge_idx][k]; 01407 int k1 = (k > 1)?(k + 1):k; 01408 tetra.vertex(k1) = h_element->vertex[k0]->index; 01409 tetra.boundary(k) = h_element->boundary[k0]->index; 01410 } 01411 tetra.vertex(2) = refined_edge->child[0]->vertex[1]->index; 01412 01413 tetra.boundaryMark() = h_element->bmark; 01414 n_twin_tetrahedron ++; 01415 } else { // 01416 Assert((n_nonactive_neighbour == 0 && n_twin_triangle_surface == 0), ExcInternalError()); 01417 01418 tetra.vertex().resize(4); 01419 tetra.boundary().resize(4); 01420 for (int k = 0; k < 4;++ k) { 01421 tetra.vertex(k) = h_element->vertex[k]->index; 01422 tetra.boundary(k) = h_element->boundary[k]->index; 01423 } 01424 01425 tetra.boundaryMark() = h_element->bmark; 01426 n_tetrahedron ++; 01427 } 01428 } 01429 01430 geometry_tree->unlock(); 01431 std::cerr << "\n\tnodes: " << n_node 01432 << "; edges: " << n_edge 01433 << "; surface: " << n_surface 01434 << "; elements: " << n_element 01435 << " (" << n_tetrahedron 01436 << ", " << n_twin_tetrahedron 01437 << ", " << n_four_tetrahedron 01438 << ")" << std::endl; 01439 01440 if (renumerate) renumerateElement(); 01441 } 01442 01443 01444 template <> 01445 void RegularMesh<2, DOW>::writeEasyMesh(const std::string& filename) const 01446 { 01447 unsigned int i, j, k, ii[] = {0, 1, 2, 0, 1, 2, 0, 1, 2}; 01448 unsigned int n_node = n_geometry(0); 01449 unsigned int n_side = n_geometry(1); 01450 unsigned int n_element = n_geometry(2); 01451 unsigned int n_triangle = 0; 01452 unsigned int n_twin_triangle = 0; 01453 std::cerr << "Writing easymesh data file ..." << std::endl; 01454 std::cerr << " preparing data ..." << std::flush; 01455 std::vector<int> index(n_element, -1); 01456 for (i = 0;i < n_element;i ++) { 01457 j = geometry(2, i).n_vertex(); 01458 if (j == 3) 01459 n_triangle ++; 01460 else if (j == 4) 01461 index[i] = n_twin_triangle ++; 01462 else 01463 Assert(false, ExcInternalError()); 01464 } 01465 01467 std::vector<std::vector<int> > side_neighbour 01468 (2, std::vector<int>(n_side, -1)); 01469 for (i = 0;i < n_element;i ++) { 01470 const GeometryBM& the_ele = geometry(2, i); 01471 if (index[i] == -1) { // this is a triangle 01472 for (j = 0;j < 3;j ++) { 01473 k = the_ele.boundary(j); 01474 const GeometryBM& the_side = geometry(1, k); 01475 if (the_side.vertex(0) == the_ele.vertex(ii[j+1])) { 01476 Assert(the_side.vertex(1) == the_ele.vertex(ii[j+2]), ExcInternalError()); 01477 side_neighbour[1][k] = i; 01478 } 01479 else if (the_side.vertex(0) == the_ele.vertex(ii[j+2])) { 01480 Assert(the_side.vertex(1) == the_ele.vertex(ii[j+1]), ExcInternalError()); 01481 side_neighbour[0][k] = i; 01482 } 01483 else { 01484 Assert(false, ExcInternalError()); // something must be wrong! 01485 } 01486 } 01487 } 01488 else { // this is a twin triangle 01489 // the 0-th side 01490 k = the_ele.boundary(0); 01491 const GeometryBM& the_side_0 = geometry(1, k); 01492 if (the_side_0.vertex(0) == the_ele.vertex(0)) { 01493 Assert(the_side_0.vertex(1) == the_ele.vertex(1), ExcInternalError()); 01494 side_neighbour[1][k] = i; 01495 } 01496 else if (the_side_0.vertex(0) == the_ele.vertex(1)) { 01497 Assert(the_side_0.vertex(1) == the_ele.vertex(0), ExcInternalError()); 01498 side_neighbour[0][k] = i; 01499 } 01500 else { 01501 Assert(false, ExcInternalError()); 01502 } 01503 // the 1-st side 01504 k = the_ele.boundary(1); 01505 const GeometryBM& the_side_1 = geometry(1, k); 01506 if (the_side_1.vertex(0) == the_ele.vertex(1)) { 01507 Assert(the_side_1.vertex(1) == the_ele.vertex(2), ExcInternalError()); 01508 side_neighbour[1][k] = i; 01509 } 01510 else if (the_side_1.vertex(0) == the_ele.vertex(2)) { 01511 Assert(the_side_1.vertex(1) == the_ele.vertex(1), ExcInternalError()); 01512 side_neighbour[0][k] = i; 01513 } 01514 else { 01515 Assert(false, ExcInternalError()); 01516 } 01517 // the 2-nd side 01518 k = the_ele.boundary(2); 01519 const GeometryBM& the_side_2 = geometry(1, k); 01520 if (the_side_2.vertex(0) == the_ele.vertex(2)) { 01521 Assert(the_side_2.vertex(1) == the_ele.vertex(3), ExcInternalError()); 01522 side_neighbour[1][k] = n_element + index[i]; 01523 } 01524 else if (the_side_2.vertex(0) == the_ele.vertex(3)) { 01525 Assert(the_side_2.vertex(1) == the_ele.vertex(2), ExcInternalError()); 01526 side_neighbour[0][k] = n_element + index[i]; 01527 } 01528 else { 01529 Assert(false, ExcInternalError()); 01530 } 01531 // the 3th side 01532 k = the_ele.boundary(3); 01533 const GeometryBM& the_side_3 = geometry(1, k); 01534 if (the_side_3.vertex(0) == the_ele.vertex(3)) { 01535 Assert(the_side_3.vertex(1) == the_ele.vertex(0), ExcInternalError()); 01536 side_neighbour[1][k] = n_element + index[i]; 01537 } 01538 else if (the_side_3.vertex(0) == the_ele.vertex(0)) { 01539 Assert(the_side_3.vertex(1) == the_ele.vertex(3), ExcInternalError()); 01540 side_neighbour[0][k] = n_element + index[i]; 01541 } 01542 else { 01543 Assert(false, ExcInternalError()); 01544 } 01545 } 01546 } 01547 01549 std::vector<std::vector<int> > element_neighbour 01550 (3, std::vector<int>(n_element + n_twin_triangle, -1)); 01551 for (i = 0;i < n_element;i ++) { 01552 const GeometryBM& the_ele = geometry(2, i); 01553 if (index[i] == -1) { 01554 for (j = 0;j < 3;j ++) { 01555 k = the_ele.boundary(j); 01556 const GeometryBM& the_side = geometry(1, k); 01557 if (the_side.vertex(0) == the_ele.vertex(ii[j+1])) { 01558 element_neighbour[j][i] = side_neighbour[0][k]; 01559 } 01560 else { // if (the_side->vertex[0] == the_ele->vertex[ii[j+2]]) 01561 element_neighbour[j][i] = side_neighbour[1][k]; 01562 } 01563 } 01564 } 01565 else { 01566 element_neighbour[1][i] = n_element + index[i]; 01567 element_neighbour[2][n_element + index[i]] = i; 01568 k = the_ele.boundary(0); 01569 const GeometryBM& the_side_0 = geometry(1, k); 01570 if (the_side_0.vertex(0) == the_ele.vertex(0)) 01571 element_neighbour[2][i] = side_neighbour[0][k]; 01572 else 01573 element_neighbour[2][i] = side_neighbour[1][k]; 01574 k = the_ele.boundary(1); 01575 const GeometryBM& the_side_1 = geometry(1, k); 01576 if (the_side_1.vertex(0) == the_ele.vertex(1)) 01577 element_neighbour[0][i] = side_neighbour[0][k]; 01578 else 01579 element_neighbour[0][i] = side_neighbour[1][k]; 01580 k = the_ele.boundary(2); 01581 const GeometryBM& the_side_2 = geometry(1, k); 01582 if (the_side_2.vertex(0) == the_ele.vertex(2)) 01583 element_neighbour[0][n_element + index[i]] = side_neighbour[0][k]; 01584 else 01585 element_neighbour[0][n_element + index[i]] = side_neighbour[1][k]; 01586 k = the_ele.boundary(3); 01587 const GeometryBM& the_side_3 = geometry(1, k); 01588 if (the_side_3.vertex(0) == the_ele.vertex(3)) 01589 element_neighbour[1][n_element + index[i]] = side_neighbour[0][k]; 01590 else 01591 element_neighbour[1][n_element + index[i]] = side_neighbour[1][k]; 01592 } 01593 } 01594 std::cerr << " OK!" << std::endl; 01595 01596 std::cerr << " writing node data ..." << std::flush; 01597 std::ofstream os((filename + ".n").c_str()); 01598 os << n_node << "\t" 01599 << n_element + n_twin_triangle << "\t" 01600 << n_side + n_twin_triangle << " **(Nnd, Nee, Nsd)**\n"; 01601 os.precision(8); 01602 os.setf(std::ios::scientific, std::ios::floatfield); 01603 for (i = 0;i < n_node;i ++) { 01604 j = geometry(0, i).vertex(0); 01605 os << i << "\t" 01606 << point(j)[0] << "\t" 01607 << point(j)[1] << "\t" 01608 << boundaryMark(0, i) << "\n"; 01609 } 01610 os << "---------------------------------------------------------------\n" 01611 << " n: x y \n"; 01612 os.close(); 01613 std::cerr << " OK!" << std::endl; 01614 01615 std::cerr << " writing side data ..." << std::flush; 01616 os.open((filename + ".s").c_str()); 01617 os << n_side + n_twin_triangle << "\n"; 01618 for (i = 0;i < n_side;i ++) { 01619 const GeometryBM& the_side = geometry(1, i); 01620 os << i << "\t" 01621 << the_side.vertex(0) << "\t" 01622 << the_side.vertex(1) << "\t" 01623 << side_neighbour[0][i] << "\t" 01624 << side_neighbour[1][i] << "\t" 01625 << boundaryMark(1, i) << "\n"; 01626 } 01627 for (i = 0;i < n_element;i ++) { 01628 if (index[i] == -1) continue; 01629 os << n_side + index[i] << "\t" 01630 << geometry(2,i).vertex(0) << "\t" 01631 << geometry(2,i).vertex(2) << "\t" 01632 << i << "\t" 01633 << n_element + index[i] << "\t" 01634 << boundaryMark(2, i) << "\n"; 01635 } 01636 os << "---------------------------------------------------------------\n" 01637 << " s: c d ea eb \n"; 01638 os.close(); 01639 std::cerr << " OK!" << std::endl; 01640 01641 std::cerr << " writing element data ..." << std::flush; 01642 os.open((filename+".e").c_str()); 01643 os << n_element + n_twin_triangle << "\t" 01644 << n_node << "\t" 01645 << n_side + n_twin_triangle << " **(Nee, Nnd, Nsd)**\n"; 01646 for (i = 0;i < n_element;i ++) { 01647 const GeometryBM& the_ele = geometry(2, i); 01648 if (index[i] == -1) { 01649 os << i << "\t" 01650 << the_ele.vertex(0) << "\t" 01651 << the_ele.vertex(1) << "\t" 01652 << the_ele.vertex(2) << "\t" 01653 << element_neighbour[0][i] << "\t" 01654 << element_neighbour[1][i] << "\t" 01655 << element_neighbour[2][i] << "\t" 01656 << the_ele.boundary(0) << "\t" 01657 << the_ele.boundary(1) << "\t" 01658 << the_ele.boundary(2) << "\n"; 01659 } 01660 else { 01661 os << i << "\t" 01662 << the_ele.vertex(0) << "\t" 01663 << the_ele.vertex(1) << "\t" 01664 << the_ele.vertex(2) << "\t" 01665 << element_neighbour[0][i] << "\t" 01666 << element_neighbour[1][i] << "\t" 01667 << element_neighbour[2][i] << "\t" 01668 << the_ele.boundary(1) << "\t" 01669 << n_side + index[i] << "\t" 01670 << the_ele.boundary(0) << "\n"; 01671 } 01672 } 01673 for (i = 0;i < n_element;i ++) { 01674 if (index[i] == -1) continue; 01675 const GeometryBM& the_ele = geometry(2, i); 01676 j = n_element + index[i]; 01677 os << j << "\t" 01678 << geometry(2,i).vertex(0) << "\t" 01679 << geometry(2,i).vertex(2) << "\t" 01680 << geometry(2,i).vertex(3) << "\t" 01681 << element_neighbour[0][j] << "\t" 01682 << element_neighbour[1][j] << "\t" 01683 << element_neighbour[2][j] << "\t" 01684 << the_ele.boundary(2) << "\t" 01685 << the_ele.boundary(3) << "\t" 01686 << n_side + index[i] << "\n"; 01687 } 01688 os << "---------------------------------------------------------------\n" 01689 << " e: i, j, k, ei, ej, ek, si, sj, sk \n"; 01690 os.close(); 01691 std::cerr << " OK!" << std::endl; 01692 } 01693 01694 01695 01696 template <> 01697 void RegularMesh<2, DOW>::writeTecplotData(const std::string& filename) const 01698 { 01699 int i; 01700 std::cerr << "Write mesh data into Tecplot data file " 01701 << filename << " ... " << std::flush; 01702 std::ofstream os(filename.c_str()); 01703 os << "TITLE = \x22" << "2D mesh data generated by AFEPack" << "\x22\n" 01704 << "VARIABLES = \x22" << "X\x22, \x22" << "Y\x22\n"; 01705 os.precision(8); 01706 os.setf(std::ios::scientific, std::ios::floatfield); 01707 01708 int n_node = n_point(); 01709 int n_element = n_geometry(2); 01710 os << "ZONE N=" << n_node 01711 << ",E=" << n_element 01712 << ",F=FEPOINT ET=TRIANGLE\n"; 01713 for (i = 0;i < n_node;i ++) 01714 os << point(i) << "\n"; 01715 os << "\n"; 01716 for (i = 0;i < n_element;i ++) { 01717 int n_vertex = geometry(2,i).n_vertex(); 01718 switch (n_vertex) { 01719 case 3: 01720 os << geometry(0,geometry(2,i).vertex(0)).vertex(0)+1 << "\t" 01721 << geometry(0,geometry(2,i).vertex(1)).vertex(0)+1 << "\t" 01722 << geometry(0,geometry(2,i).vertex(2)).vertex(0)+1 << "\n"; 01723 break; 01724 case 4: 01725 os << geometry(0,geometry(2,i).vertex(0)).vertex(0)+1 << "\t" 01726 << geometry(0,geometry(2,i).vertex(1)).vertex(0)+1 << "\t" 01727 << geometry(0,geometry(2,i).vertex(2)).vertex(0)+1 << "\n"; 01728 break; 01729 default: 01730 Assert(false, ExcInternalError()); 01731 } 01732 } 01733 os.close(); 01734 std::cerr << "OK!" << std::endl; 01735 } 01736 01737 01738 01739 template <> 01740 void RegularMesh<3, DOW>::writeTecplotData(const std::string& filename) const 01741 { 01742 int i; 01743 std::cerr << "Write mesh data into Tecplot data file " 01744 << filename << " ... " << std::flush; 01745 std::ofstream os(filename.c_str()); 01746 os << "TITLE = \x22" << "3D mesh data generated by AFEPack" << "\x22\n" 01747 << "VARIABLES = \x22" << "X\x22, \x22" << "Y\x22, \x22" << "Z\x22\n"; 01748 os.precision(8); 01749 os.setf(std::ios::scientific, std::ios::floatfield); 01750 01751 int n_node = n_point(); 01752 int n_element = n_geometry(3); 01753 os << "ZONE N=" << n_node 01754 << ",E=" << n_element 01755 << ",F=FEPOINT ET=TETRAHEDRON\n"; 01756 for (i = 0;i < n_node;i ++) 01757 os << point(i) << "\n"; 01758 os << "\n"; 01759 for (i = 0;i < n_element;i ++) { 01760 int n_vertex = geometry(3,i).n_vertex(); 01761 switch (n_vertex) { 01762 case 4: 01763 os << geometry(0,geometry(3,i).vertex(0)).vertex(0) + 1 << "\t" 01764 << geometry(0,geometry(3,i).vertex(1)).vertex(0) + 1 << "\t" 01765 << geometry(0,geometry(3,i).vertex(2)).vertex(0) + 1 << "\t" 01766 << geometry(0,geometry(3,i).vertex(3)).vertex(0) + 1 << "\n"; 01767 break; 01768 case 5: 01769 os << geometry(0,geometry(3,i).vertex(0)).vertex(0) + 1 << "\t" 01770 << geometry(0,geometry(3,i).vertex(1)).vertex(0) + 1 << "\t" 01771 << geometry(0,geometry(3,i).vertex(3)).vertex(0) + 1 << "\t" 01772 << geometry(0,geometry(3,i).vertex(4)).vertex(0) + 1 << "\n"; 01773 break; 01774 case 7: 01775 os << geometry(0,geometry(3,i).vertex(0)).vertex(0) + 1 << "\t" 01776 << geometry(0,geometry(3,i).vertex(1)).vertex(0) + 1 << "\t" 01777 << geometry(0,geometry(3,i).vertex(2)).vertex(0) + 1 << "\t" 01778 << geometry(0,geometry(3,i).vertex(3)).vertex(0) + 1 << "\n"; 01779 break; 01780 default: 01781 Assert(false, ExcInternalError()); 01782 } 01783 } 01784 os.close(); 01785 std::cerr << "OK!" << std::endl; 01786 } 01787 01788 template <> 01789 void RegularMesh<1, DOW>::writeOpenDXData(const std::string& filename) const 01790 { 01791 std::cout << "Not implemented." << std::endl; 01792 } 01793 01794 template <> 01795 void RegularMesh<2, DOW>::writeOpenDXData(const std::string& filename) const 01796 { 01797 std::ofstream os(filename.c_str()); 01798 01799 os.precision(8); 01800 os.setf(std::ios::scientific, std::ios::floatfield); 01801 int n_node = n_point(); 01802 int i, j; 01803 01804 os << "object 1 class array type float rank 1 shape " << DOW << " item " 01805 << n_node << " data follows\n"; 01806 for (i = 0;i < n_node;i ++) { 01807 os << point(geometry(0,i).vertex(0)) << "\n"; 01808 } 01809 01810 int n_element = n_geometry(2); 01811 for (i = 0, j = 0;i < n_element;i ++) { 01812 switch (geometry(2,i).n_vertex()) { 01813 case 3: 01814 j += 1; 01815 break; 01816 case 4: 01817 j += 2; 01818 break; 01819 default: 01820 Assert(false, ExcInternalError()); 01821 } 01822 } 01823 os << "\nobject 2 class array type int rank 1 shape 3 item " 01824 << j << " data follows\n"; 01825 for (i = 0;i < n_element;i ++) { 01826 switch (geometry(2,i).n_vertex()) { 01827 case 3: 01828 os << geometry(2,i).vertex(0) << "\t" 01829 << geometry(2,i).vertex(1) << "\t" 01830 << geometry(2,i).vertex(2) << "\t\n"; 01831 break; 01832 case 4: 01833 os << geometry(2,i).vertex(0) << "\t" 01834 << geometry(2,i).vertex(1) << "\t" 01835 << geometry(2,i).vertex(2) << "\t\n"; 01836 os << geometry(2,i).vertex(0) << "\t" 01837 << geometry(2,i).vertex(2) << "\t" 01838 << geometry(2,i).vertex(3) << "\t\n"; 01839 break; 01840 default: 01841 Assert(false, ExcInternalError()); 01842 } 01843 } 01844 os << "attribute \"element type\" string \"triangles\"\n" 01845 << "attribute \"ref\" string \"positions\"\n\n"; 01846 01847 os << "object \"FEMFunction-2d\" class field\n" 01848 << "component \"positions\" value 1\n" 01849 << "component \"connections\" value 2\n" 01850 << "end\n"; 01851 os.close(); 01852 } 01853 01854 01855 01856 template <> 01857 void RegularMesh<3, DOW>::writeOpenDXData(const std::string& filename) const 01858 { 01859 std::ofstream os(filename.c_str()); 01860 01861 os.precision(8); 01862 os.setf(std::ios::scientific, std::ios::floatfield); 01863 int n_node = n_point(); 01864 int i, j; 01865 01866 os << "object 1 class array type float rank 1 shape " << DOW << " item " 01867 << n_node << " data follows\n"; 01868 for (i = 0;i < n_node;i ++) { 01869 os << point(geometry(0,i).vertex(0)) << "\n"; 01870 } 01871 01872 int n_element = n_geometry(3); 01873 for (i = 0, j = 0;i < n_element;i ++) { 01874 switch (geometry(3,i).n_vertex()) { 01875 case 4: 01876 j += 1; 01877 break; 01878 case 5: 01879 j += 2; 01880 break; 01881 case 7: 01882 j += 4; 01883 break; 01884 default: 01885 Assert(false, ExcInternalError()); 01886 } 01887 } 01888 os << "\nobject 2 class array type int rank 1 shape 4 item " 01889 << j << " data follows\n"; 01890 for (i = 0;i < n_element;i ++) { 01891 switch (geometry(3,i).n_vertex()) { 01892 case 4: 01893 os << geometry(3,i).vertex(0) << "\t" 01894 << geometry(3,i).vertex(1) << "\t" 01895 << geometry(3,i).vertex(2) << "\t" 01896 << geometry(3,i).vertex(3) << "\t\n"; 01897 break; 01898 case 5: 01899 os << geometry(3,i).vertex(0) << "\t" 01900 << geometry(3,i).vertex(1) << "\t" 01901 << geometry(3,i).vertex(2) << "\t" 01902 << geometry(3,i).vertex(4) << "\t\n"; 01903 os << geometry(3,i).vertex(0) << "\t" 01904 << geometry(3,i).vertex(2) << "\t" 01905 << geometry(3,i).vertex(3) << "\t" 01906 << geometry(3,i).vertex(4) << "\t\n"; 01907 break; 01908 case 7: 01909 os << geometry(3,i).vertex(0) << "\t" 01910 << geometry(3,i).vertex(1) << "\t" 01911 << geometry(3,i).vertex(6) << "\t" 01912 << geometry(3,i).vertex(5) << "\t\n"; 01913 os << geometry(3,i).vertex(0) << "\t" 01914 << geometry(3,i).vertex(2) << "\t" 01915 << geometry(3,i).vertex(4) << "\t" 01916 << geometry(3,i).vertex(6) << "\t\n"; 01917 os << geometry(3,i).vertex(0) << "\t" 01918 << geometry(3,i).vertex(3) << "\t" 01919 << geometry(3,i).vertex(5) << "\t" 01920 << geometry(3,i).vertex(4) << "\t\n"; 01921 os << geometry(3,i).vertex(0) << "\t" 01922 << geometry(3,i).vertex(4) << "\t" 01923 << geometry(3,i).vertex(5) << "\t" 01924 << geometry(3,i).vertex(6) << "\t\n"; 01925 break; 01926 default: 01927 Assert(false, ExcInternalError()); 01928 } 01929 } 01930 os << "attribute \"element type\" string \"tetrahedra\"\n" 01931 << "attribute \"ref\" string \"positions\"\n\n"; 01932 01933 os << "object \"FEMFunction-3d\" class field\n" 01934 << "component \"positions\" value 1\n" 01935 << "component \"connections\" value 2\n" 01936 << "end\n"; 01937 os.close(); 01938 } 01939 01940 01941 01942 template <> 01943 void HGeometryTree<2, DOW>::readEasyMesh(const std::string& filename) 01944 { 01945 std::cerr << "Reading easymesh data file ..." << std::endl; 01946 01947 std::ifstream is((filename + ".n").c_str()); 01948 int i, j, k, l; 01949 int n_node, n_side, n_element; 01950 char dummy[64]; 01951 is >> n_node >> n_element >> n_side; 01952 is.getline(dummy, 64); 01953 std::vector<HGeometry<0, DOW> *> node(n_node); 01954 std::vector<HGeometry<1, DOW> *> side(n_side); 01955 std::vector<HGeometry<2, DOW> *> element(n_element); 01956 std::cerr << "\treading the nodes data ..." << std::flush; 01957 for (i = 0;i < n_node;i ++) { 01958 node[i] = new HGeometry<0, DOW>(); 01959 Assert(node[i] != NULL, ExcOutOfMemory()); 01960 is >> j 01961 >> *(dynamic_cast<afepack::Point<DOW> *>(node[i])) 01962 >> node[i]->bmark; 01963 } 01964 is.close(); 01965 std::cerr << " OK!" << std::endl; 01966 01967 is.open((filename + ".s").c_str()); 01968 is >> i; 01969 Assert(i == n_side, ExcInternalError()); 01970 std::cerr << "\treading the sides data ..." << std::flush; 01971 for (i = 0;i < n_side;i ++) { 01972 side[i] = new HGeometry<1, DOW>(); 01973 Assert(side[i] != NULL, ExcOutOfMemory()); 01974 is >> l >> j >> k; 01975 side[i]->vertex[0] = node[j]; 01976 side[i]->vertex[1] = node[k]; 01977 is >> j >> k >> side[i]->bmark; 01978 } 01979 is.close(); 01980 std::cerr << " OK!" << std::endl; 01981 01982 is.open((filename + ".e").c_str()); 01983 is >> i >> j >> k; 01984 Assert(i == n_element, ExcInternalError()); 01985 Assert(j == n_node, ExcInternalError()); 01986 Assert(k == n_side, ExcInternalError()); 01987 is.getline(dummy, 64); 01988 std::cerr << "\treading the elements data ..." << std::flush; 01989 for (i = 0;i < n_element;i ++) { 01990 element[i] = new HGeometry<2, DOW>(); 01991 Assert(element[i] != NULL, ExcOutOfMemory()); 01992 is >> l >> j >> k >> l; 01993 element[i]->vertex[0] = node[j]; 01994 element[i]->vertex[1] = node[k]; 01995 element[i]->vertex[2] = node[l]; 01996 element[i]->bmark = 0; 01997 is >> j >> k >> l; 01998 is >> j >> k >> l; 01999 element[i]->boundary[0] = side[j]; 02000 element[i]->boundary[1] = side[k]; 02001 element[i]->boundary[2] = side[l]; 02002 root_element.push_back(element[i]); 02003 } 02004 is.close(); 02005 std::cerr << " OK!" << std::endl; 02006 } 02007 02008 02009 template <> 02010 void RegularMesh<1, DOW>::writeSimplestSimplexMesh(const std::string& filename) const 02011 { 02012 std::cout << "Not implemented." << std::endl; 02013 } 02014 02015 template <> 02016 void RegularMesh<2, DOW>::writeSimplestSimplexMesh(const std::string& filename) const 02017 { 02018 std::ofstream os(filename.c_str()); 02019 02020 os.precision(8); 02021 os.setf(std::ios::scientific, std::ios::floatfield); 02022 int n_node = n_point(); 02023 int i, j; 02024 02025 os << n_node << "\n"; 02026 for (i = 0;i < n_node;i ++) { 02027 os << point(geometry(0,i).vertex(0)) << "\n"; 02028 } 02029 02030 int n_element = n_geometry(2); 02031 for (i = 0, j = 0;i < n_element;i ++) { 02032 switch (geometry(2,i).n_vertex()) { 02033 case 3: 02034 j += 1; 02035 break; 02036 case 4: 02037 j += 2; 02038 break; 02039 default: 02040 Assert(false, ExcInternalError()); 02041 } 02042 } 02043 os << j << "\n"; 02044 for (i = 0;i < n_element;i ++) { 02045 switch (geometry(2,i).n_vertex()) { 02046 case 3: 02047 os << geometry(2,i).vertex(0) << "\t" 02048 << geometry(2,i).vertex(1) << "\t" 02049 << geometry(2,i).vertex(2) << "\t\n"; 02050 break; 02051 case 4: 02052 os << geometry(2,i).vertex(0) << "\t" 02053 << geometry(2,i).vertex(1) << "\t" 02054 << geometry(2,i).vertex(2) << "\t\n"; 02055 os << geometry(2,i).vertex(0) << "\t" 02056 << geometry(2,i).vertex(2) << "\t" 02057 << geometry(2,i).vertex(3) << "\t\n"; 02058 break; 02059 default: 02060 Assert(false, ExcInternalError()); 02061 } 02062 } 02063 os.close(); 02064 } 02065 02066 02067 02068 template <> 02069 void RegularMesh<3, DOW>::writeSimplestSimplexMesh(const std::string& filename) const 02070 { 02071 std::ofstream os(filename.c_str()); 02072 02073 os.precision(8); 02074 os.setf(std::ios::scientific, std::ios::floatfield); 02075 int n_node = n_point(); 02076 int i, j; 02077 02078 os << n_node << "\n"; 02079 for (i = 0;i < n_node;i ++) { 02080 os << point(geometry(0,i).vertex(0)) << "\n"; 02081 } 02082 02083 int n_element = n_geometry(3); 02084 for (i = 0, j = 0;i < n_element;i ++) { 02085 switch (geometry(3,i).n_vertex()) { 02086 case 4: 02087 j += 1; 02088 break; 02089 case 5: 02090 j += 2; 02091 break; 02092 case 7: 02093 j += 4; 02094 break; 02095 default: 02096 Assert(false, ExcInternalError()); 02097 } 02098 } 02099 os << j << "\n"; 02100 for (i = 0;i < n_element;i ++) { 02101 switch (geometry(3,i).n_vertex()) { 02102 case 4: 02103 os << geometry(3,i).vertex(0) << "\t" 02104 << geometry(3,i).vertex(1) << "\t" 02105 << geometry(3,i).vertex(2) << "\t" 02106 << geometry(3,i).vertex(3) << "\t\n"; 02107 break; 02108 case 5: 02109 os << geometry(3,i).vertex(0) << "\t" 02110 << geometry(3,i).vertex(1) << "\t" 02111 << geometry(3,i).vertex(2) << "\t" 02112 << geometry(3,i).vertex(4) << "\t\n"; 02113 os << geometry(3,i).vertex(0) << "\t" 02114 << geometry(3,i).vertex(2) << "\t" 02115 << geometry(3,i).vertex(3) << "\t" 02116 << geometry(3,i).vertex(4) << "\t\n"; 02117 break; 02118 case 7: 02119 os << geometry(3,i).vertex(0) << "\t" 02120 << geometry(3,i).vertex(1) << "\t" 02121 << geometry(3,i).vertex(6) << "\t" 02122 << geometry(3,i).vertex(5) << "\t\n"; 02123 os << geometry(3,i).vertex(0) << "\t" 02124 << geometry(3,i).vertex(2) << "\t" 02125 << geometry(3,i).vertex(4) << "\t" 02126 << geometry(3,i).vertex(6) << "\t\n"; 02127 os << geometry(3,i).vertex(0) << "\t" 02128 << geometry(3,i).vertex(3) << "\t" 02129 << geometry(3,i).vertex(5) << "\t" 02130 << geometry(3,i).vertex(4) << "\t\n"; 02131 os << geometry(3,i).vertex(0) << "\t" 02132 << geometry(3,i).vertex(4) << "\t" 02133 << geometry(3,i).vertex(5) << "\t" 02134 << geometry(3,i).vertex(6) << "\t\n"; 02135 break; 02136 default: 02137 Assert(false, ExcInternalError()); 02138 } 02139 } 02140 os.close(); 02141 } 02142 02143 template <> 02144 void RegularMesh<1, DOW>::writeSimplexMesh(const std::string& filename) const 02145 { 02146 std::cout << "Not implemented." << std::endl; 02147 } 02148 02149 template <> 02150 void RegularMesh<2, DOW>::writeSimplexMesh(const std::string& filename) const 02151 { 02152 std::ofstream os(filename.c_str()); 02153 02154 os.precision(8); 02155 os.setf(std::ios::scientific, std::ios::floatfield); 02156 int n_node = n_point(); 02157 int i, j; 02158 02159 os << n_node << "\n"; 02160 for (i = 0;i < n_node;i ++) { 02161 os << point(geometry(0,i).vertex(0)) << "\n"; 02162 } 02163 02164 os << "\n" << n_node << "\n"; 02165 for (i = 0;i < n_node;i ++) { 02166 os << i << "\n" 02167 << "1\t" << i << "\n" 02168 << "1\t" << i << "\n" 02169 << geometry(0,i).boundaryMark() << "\n"; 02170 } 02171 02172 int n_twin_triangle = 0; 02173 int n_element = n_geometry(2); 02174 for (i = 0;i < n_element;i ++) { 02175 if (geometry(2,i).n_vertex() == 4) { 02176 n_twin_triangle += 1; 02177 } 02178 } 02179 02180 int n_side = n_geometry(1); 02181 os << "\n" << n_side + n_twin_triangle << "\n"; 02182 for (i = 0;i < n_side;i ++) { 02183 os << i << "\n" 02184 << "2\t" << geometry(1,i).vertex(0) << " " << geometry(1,i).vertex(1) << "\n" 02185 << "2\t" << geometry(1,i).vertex(0) << " " << geometry(1,i).vertex(1) << "\n" 02186 << geometry(1,i).boundaryMark() << "\n"; 02187 } 02188 for (i = 0, j = 0;i < n_element;i ++) { 02189 if (geometry(2,i).n_vertex() == 4) { 02190 os << n_side + j << "\n" 02191 << "2\t" << geometry(2,i).vertex(0) << " " << geometry(2,i).vertex(2) << "\n" 02192 << "2\t" << geometry(2,i).vertex(0) << " " << geometry(2,i).vertex(2) << "\n" 02193 << geometry(2,i).boundaryMark() << "\n"; 02194 j += 1; 02195 } 02196 } 02197 02198 os << "n" << n_element + n_twin_triangle << "\n"; 02199 for (i = 0, j = 0;i < n_element;i ++) { 02200 const GeometryBM& geo = geometry(2,i); 02201 switch (geo.n_vertex()) { 02202 case 3: 02203 os << i + j << "\n" 02204 << "3\t" << geo.vertex(0) << " " << geo.vertex(1) << " " << geo.vertex(2) << "\n" 02205 << "3\t" << geo.boundary(0) << " " << geo.boundary(1) << " " << geo.boundary(2) << "\n" 02206 << geo.boundaryMark() << "\n"; 02207 break; 02208 case 4: 02209 os << i + j << "\n" 02210 << "3\t" << geo.vertex(0) << " " << geo.vertex(1) << " " << geo.vertex(2) << "\n" 02211 << "3\t" << geo.boundary(1) << " " << n_side + j << " " << geo.boundary(0) << "\n" 02212 << geo.boundaryMark() << "\n"; 02213 j += 1; 02214 os << i + j << "\n" 02215 << "3\t" << geo.vertex(0) << " " << geo.vertex(2) << " " << geo.vertex(3) << "\n" 02216 << "3\t" << geo.boundary(2) << " " << geo.boundary(3) << " " << n_side + j << "\n" 02217 << geo.boundaryMark() << "\n"; 02218 break; 02219 default: 02220 Assert(false, ExcInternalError()); 02221 } 02222 } 02223 os.close(); 02224 } 02225 02226 02227 02228 template <> 02229 void RegularMesh<3, DOW>::writeSimplexMesh(const std::string& filename) const 02230 { 02231 std::ofstream os(filename.c_str()); 02232 02233 os.precision(8); 02234 os.setf(std::ios::scientific, std::ios::floatfield); 02235 int n_node = n_point(); 02236 int i, j; 02237 02238 os << n_node << "\n"; 02239 for (i = 0;i < n_node;i ++) { 02240 os << point(geometry(0,i).vertex(0)) << "\n"; 02241 } 02242 02243 os << "\n" << n_node << "\n"; 02244 for (i = 0;i < n_node;i ++) { 02245 os << i << "\n" 02246 << "1\t" << i << "\n" 02247 << "1\t" << i << "\n" 02248 << geometry(0,i).boundaryMark() << "\n"; 02249 } 02250 02251 int n_twin_triangle = 0; 02252 int n_side = n_geometry(1); 02253 int n_face = n_geometry(2); 02254 for (i = 0;i < n_face;++ i) { 02255 if (geometry(2,i).n_vertex() == 4) { 02256 n_twin_triangle += 1; 02257 } 02258 } 02259 02260 int n_twin_tetrahedron = 0; 02261 int n_four_tetrahedron = 0; 02262 int n_element = n_geometry(3); 02263 for (i = 0, j = 0;i < n_element;i ++) { 02264 switch (geometry(3,i).n_vertex()) { 02265 case 5: 02266 n_twin_tetrahedron += 1; 02267 break; 02268 case 7: 02269 n_four_tetrahedron += 1; 02270 break; 02271 } 02272 } 02273 02274 os << "\n" << n_side + n_twin_triangle << "\n"; 02275 for (i = 0;i < n_side;i ++) { 02276 os << i << "\n" 02277 << "2\t" << geometry(1,i).vertex(0) << " " << geometry(1,i).vertex(1) << "\n" 02278 << "2\t" << geometry(1,i).vertex(0) << " " << geometry(1,i).vertex(1) << "\n" 02279 << geometry(1,i).boundaryMark() << "\n"; 02280 } 02281 for (i = 0, j = 0;i < n_face;i ++) { 02282 if (geometry(2,i).n_vertex() == 4) { 02283 os << n_side + j << "\n" 02284 << "2\t" << geometry(2,i).vertex(0) << " " << geometry(2,i).vertex(2) << "\n" 02285 << "2\t" << geometry(2,i).vertex(0) << " " << geometry(2,i).vertex(2) << "\n" 02286 << geometry(2,i).boundaryMark() << "\n"; 02287 j += 1; 02288 } 02289 } 02290 02291 os << "n" << (n_face + 02292 n_twin_triangle + 02293 n_twin_tetrahedron + 02294 3*n_four_tetrahedron) << "\n"; 02295 for (i = 0;i < n_face;i ++) { 02296 const GeometryBM& geo = geometry(2,i); 02297 switch (geo.n_vertex()) { 02298 case 3: 02299 os << i + j << "\n" 02300 << "3\t" << geo.vertex(0) << " " << geo.vertex(1) << " " << geo.vertex(2) << "\n" 02301 << "3\t" << geo.boundary(0) << " " << geo.boundary(1) << " " << geo.boundary(2) << "\n" 02302 << geo.boundaryMark() << "\n"; 02303 break; 02304 case 4: 02305 os << i + j << "\n" 02306 << "3\t" << geo.vertex(0) << " " << geo.vertex(1) << " " << geo.vertex(2) << "\n" 02307 << "3\t" << geo.boundary(1) << " " << n_side + j << " " << geo.boundary(0) << "\n" 02308 << geo.boundaryMark() << "\n"; 02309 break; 02310 default: 02311 Assert(false, ExcInternalError()); 02312 } 02313 } 02314 for (i = 0, j = 0;i < n_face;i ++) { 02315 const GeometryBM& geo = geometry(2,i); 02316 if (geo.n_vertex() == 3) continue; 02317 os << n_face + j << "\n" 02318 << "3\t" << geo.vertex(0) << " " << geo.vertex(2) << " " << geo.vertex(3) << "\n" 02319 << "3\t" << geo.boundary(2) << " " << geo.boundary(3) << " " << n_side + j << "\n" 02320 << geo.boundaryMark() << "\n"; 02321 j += 1; 02322 } 02323 for (i = 0, j = 0;i < n_element;i ++) { 02324 const GeometryBM& geo = geometry(3,i); 02325 if (geo.n_vertex() != 5) continue; 02326 os << n_face + n_twin_triangle + j << "\n" 02327 << "3\t" << geo.vertex(0) << " " << geo.vertex(2) << " " << geo.vertex(4) << "\n" 02328 << "3\t" << n_side + geo.boundary(0) << " " 02329 //<< geo.??? 02330 << n_side + geo.boundary(3) << "\n" 02331 << geo.boundaryMark() << "\n"; 02332 j += 1; 02333 } 02334 02336 02337 os.close(); 02338 } 02339