AFEPack
HGeometry.templates.nd.h
浏览该文件的文档。
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] = &mp;
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