#include #include "meshing.hpp" namespace netgen { #include "hpref_segm.hpp" #include "hpref_trig.hpp" #include "hpref_quad.hpp" #include "hpref_tet.hpp" #include "hpref_prism.hpp" #include "hpref_hex.hpp" #include "hpref_pyramid.hpp" #include "classifyhpel.hpp" void HPRefElement :: Reset(void) { np = 8; for (int i = 0; i < 8; i++) { pnums[i] = -1; param[i][0] = param[i][1] = param[i][2] = 0; } domin=-1; domout=-1; // he: levelx = 0; levely = 0; levelz = 0; } HPRefElement :: HPRefElement () { Reset(); } HPRefElement :: HPRefElement(Element & el) : type(HP_NONE), index(el.GetIndex()), levelx(0), levely(0), levelz(0), np(el.GetNV()), domin(-1), domout(-1) //domin,out for segments { //Reset(); for (int i=0; i & GetHPRegistry()\ { static std::map registry; return registry; } HPRef_Struct * Get_HPRef_Struct (HPREF_ELEMENT_TYPE type) { HPRef_Struct * hps = NULL; /* if (type >= 100 && type < 1000) if (type != HP_TET_1F_0E_1VA) type = HP_NONETET; */ switch (type) { case HP_SEGM: hps = &refsegm; break; case HP_SEGM_SINGCORNERL: hps = &refsegm_scl; break; case HP_SEGM_SINGCORNERR: hps = &refsegm_scr; break; case HP_SEGM_SINGCORNERS: hps = &refsegm_sc2; break; case HP_TRIG: hps = &reftrig; break; case HP_TRIG_SINGCORNER: hps = &reftrig_singcorner; break; case HP_TRIG_SINGCORNER12: hps = &reftrig_singcorner12; break; case HP_TRIG_SINGCORNER123: hps = &reftrig_singcorner123; break; case HP_TRIG_SINGCORNER123_2D: hps = &reftrig_singcorner123_2D; break; case HP_TRIG_SINGEDGE: hps = &reftrig_singedge; break; case HP_TRIG_SINGEDGECORNER1: hps = &reftrig_singedgecorner1; break; case HP_TRIG_SINGEDGECORNER2: hps = &reftrig_singedgecorner2; break; case HP_TRIG_SINGEDGECORNER12: hps = &reftrig_singedgecorner12; break; case HP_TRIG_SINGEDGECORNER3: hps = &reftrig_singedgecorner3; break; case HP_TRIG_SINGEDGECORNER13: hps = &reftrig_singedgecorner13; break; case HP_TRIG_SINGEDGECORNER23: hps = &reftrig_singedgecorner23; break; case HP_TRIG_SINGEDGECORNER123: hps = &reftrig_singedgecorner123; break; case HP_TRIG_SINGEDGES: hps = &reftrig_singedges; break; case HP_TRIG_SINGEDGES2: hps = &reftrig_singedges2; break; case HP_TRIG_SINGEDGES3: hps = &reftrig_singedges3; break; case HP_TRIG_SINGEDGES23: hps = &reftrig_singedges23; break; case HP_TRIG_3SINGEDGES: hps = &reftrig_3singedges; break; case HP_TRIG_ALFELD: hps = &reftrig_Alfeld; break; case HP_TRIG_POWELL: hps = &reftrig_Powell; break; case HP_QUAD: hps = &refquad; break; case HP_DUMMY_QUAD_SINGCORNER: hps = &refdummyquad_singcorner; break; case HP_QUAD_SINGCORNER: hps = &refquad_singcorner; break; case HP_QUAD_SINGEDGE: hps = &refquad_singedge; break; case HP_QUAD_0E_2VA: hps = &refquad_0e_2va; break; case HP_QUAD_0E_2VB: hps = &refquad_0e_2vb; break; case HP_QUAD_0E_3V: hps = &refquad_0e_3v; break; case HP_QUAD_0E_4V: hps = &refquad_0e_4v; break; case HP_QUAD_1E_1VA: hps = &refquad_1e_1va; break; case HP_QUAD_1E_1VB: hps = &refquad_1e_1vb; break; case HP_QUAD_1E_1VC: hps = &refquad_1e_1vc; break; case HP_QUAD_1E_1VD: hps = &refquad_1e_1vd; break; case HP_QUAD_1E_2VA: hps = &refquad_1e_2va; break; case HP_QUAD_1E_2VB: hps = &refquad_1e_2vb; break; case HP_QUAD_1E_2VC: hps = &refquad_1e_2vc; break; case HP_QUAD_1E_2VD: hps = &refquad_1e_2vd; break; case HP_QUAD_1E_2VE: hps = &refquad_1e_2ve; break; case HP_QUAD_1E_2VF: hps = &refquad_1e_2vf; break; case HP_QUAD_1E_3VA: hps = &refquad_1e_3va; break; case HP_QUAD_1E_3VB: hps = &refquad_1e_3vb; break; case HP_QUAD_1E_3VC: hps = &refquad_1e_3vc; break; case HP_QUAD_1E_3VD: hps = &refquad_1e_3vd; break; case HP_QUAD_1E_4V: hps = &refquad_1e_4v; break; case HP_QUAD_2E: hps = &refquad_2e; break; case HP_QUAD_2E_1VA: hps = &refquad_2e_1va; break; case HP_QUAD_2E_1VB: hps = &refquad_2e_1vb; break; case HP_QUAD_2E_1VC: hps = &refquad_2e_1vc; break; case HP_QUAD_2E_2VA: hps = &refquad_2e_2va; break; case HP_QUAD_2E_2VB: hps = &refquad_2e_2vb; break; case HP_QUAD_2E_2VC: hps = &refquad_2e_2vc; break; case HP_QUAD_2E_3V: hps = &refquad_2e_3v; break; case HP_QUAD_2EB_0V: hps = &refquad_2eb_0v; break; case HP_QUAD_2EB_1VA: hps = &refquad_2eb_1va; break; case HP_QUAD_2EB_1VB: hps = &refquad_2eb_1vb; break; case HP_QUAD_2EB_2VA: hps = &refquad_2eb_2va; break; case HP_QUAD_2EB_2VB: hps = &refquad_2eb_2vb; break; case HP_QUAD_2EB_2VC: hps = &refquad_2eb_2vc; break; case HP_QUAD_2EB_2VD: hps = &refquad_2eb_2vd; break; case HP_QUAD_2EB_3VA: hps = &refquad_2eb_3va; break; case HP_QUAD_2EB_3VB: hps = &refquad_2eb_3vb; break; case HP_QUAD_2EB_4V: hps = &refquad_2eb_4v; break; case HP_QUAD_3E: hps = &refquad_3e; break; case HP_QUAD_3E_3VA: hps = &refquad_3e_3va; break; case HP_QUAD_3E_3VB: hps = &refquad_3e_3vb; break; case HP_QUAD_3E_4V: hps = &refquad_3e_4v; break; case HP_QUAD_4E: hps = &refquad_4e; break; case HP_TET: hps = &reftet; break; case HP_NONETET: hps = &refnonetet; break; case HP_TET_0E_1V: hps = &reftet_0e_1v; break; case HP_TET_0E_2V: hps = &reftet_0e_2v; break; case HP_TET_0E_3V: hps = &reftet_0e_3v; break; case HP_TET_0E_4V: hps = &reftet_0e_4v; break; case HP_TET_1E_0V: hps = &reftet_1e_0v; break; case HP_TET_1E_1VA: hps = &reftet_1e_1va; break; case HP_TET_1E_1VB: hps = &reftet_1e_1vb; break; case HP_TET_1E_2VA: hps = &reftet_1e_2va; break; case HP_TET_1E_2VB: hps = &reftet_1e_2vb; break; case HP_TET_1E_2VC: hps = &reftet_1e_2vc; break; case HP_TET_1E_2VD: hps = &reftet_1e_2vd; break; case HP_TET_1E_3VA: hps = &reftet_1e_3va; break; case HP_TET_1E_3VB: hps = &reftet_1e_3vb; break; case HP_TET_1E_4V: hps = &reftet_1e_4v; break; case HP_TET_2EA_0V: hps = &reftet_2ea_0v; break; case HP_TET_2EA_1VB: hps = &reftet_2ea_1vb; break; case HP_TET_2EA_1VC: hps = &reftet_2ea_1vc; break; case HP_TET_2EA_1VA: hps = &reftet_2ea_1va; break; case HP_TET_2EA_2VA: hps = &reftet_2ea_2va; break; case HP_TET_2EA_2VB: hps = &reftet_2ea_2vb; break; case HP_TET_2EA_2VC: hps = &reftet_2ea_2vc; break; case HP_TET_2EA_3V: hps = &reftet_2ea_3v; break; case HP_TET_2EB_0V: hps = &reftet_2eb_0v; break; case HP_TET_2EB_1V: hps = &reftet_2eb_1v; break; case HP_TET_2EB_2VA: hps = &reftet_2eb_2va; break; case HP_TET_2EB_2VB: hps = &reftet_2eb_2vb; break; case HP_TET_2EB_2VC: hps = &reftet_2eb_2vc; break; case HP_TET_2EB_3V: hps = &reftet_2eb_3v; break; case HP_TET_2EB_4V: hps = &reftet_2eb_4v; break; case HP_TET_3EA_0V: hps = &reftet_3ea_0v; break; case HP_TET_3EA_1V: hps = &reftet_3ea_1v; break; case HP_TET_3EA_2V: hps = &reftet_3ea_2v; break; case HP_TET_3EA_3V: hps = &reftet_3ea_3v; break; case HP_TET_3EB_0V: hps = &reftet_3eb_0v; break; case HP_TET_3EB_1V: hps = &reftet_3eb_1v; break; case HP_TET_3EB_2V: hps = &reftet_3eb_2v; break; case HP_TET_3EC_0V: hps = &reftet_3ec_0v; break; case HP_TET_3EC_1V: hps = &reftet_3ec_1v; break; case HP_TET_3EC_2V: hps = &reftet_3ec_2v; break; case HP_TET_3ED_3V: hps = &reftet_3ed_3v; break; case HP_TET_1F_0E_0V: hps = &reftet_1f_0e_0v; break; // case HP_TET_1F_0E_1VA: // hps = &reftet_1f_0e_1va; break; case HP_TET_1F_0E_1VB: hps = &reftet_1f_0e_1vb; break; case HP_TET_1F_0E_2V: hps = &reftet_1f_0e_2v; break; // case HP_TET_1F_0E_3V: // hps = &reftet_1f_0e_3v; break; case HP_TET_1F_1EA_0V: hps = &reftet_1f_1ea_0v; break; case HP_TET_1F_1EB_0V: hps = &reftet_1f_1eb_0v; break; case HP_TET_1F_1E_1VA: hps = &reftet_1f_1e_1va; break; case HP_TET_1F_1E_1VB: hps = &reftet_1f_1e_1vb; break; case HP_TET_1F_1E_2VA: hps = &reftet_1f_1e_2va; break; case HP_TET_1F_1E_2VB: hps = &reftet_1f_1e_2vb; break; case HP_TET_1F_2E_0VA: hps = &reftet_1f_2e_0va; break; case HP_TET_1F_2E_0VB: hps = &reftet_1f_2e_0vb; break; case HP_TET_2F_0E_0V: hps = &reftet_2f_0e_0v; break; case HP_TET_2F_0E_1V: hps = &reftet_2f_0e_1v; break; case HP_TET_2F_1E_0VA: hps = &reftet_2f_1e_0va; break; case HP_TET_2F_1E_0VB: hps = &reftet_2f_1e_0vb; break; case HP_TET_3F_0E_0V: hps = &reftet_3f_0e_0v; break; case HP_PRISM: hps = &refprism; break; case HP_PRISM_SINGEDGE: hps = &refprism_singedge; break; // case HP_PRISM_SINGEDGE_H1: // hps = &refprism_singedge_h1; break; // case HP_PRISM_SINGEDGE_H12: // hps = &refprism_singedge_h12; break; case HP_PRISM_SINGEDGE_V12: hps = &refprism_singedge_v12; break; case HP_PRISM_1FA_0E_0V: hps = &refprism_1fa_0e_0v; break; case HP_PRISM_2FA_0E_0V: hps = &refprism_2fa_0e_0v; break; case HP_PRISM_1FB_0E_0V: hps = &refprism_1fb_0e_0v; break; case HP_PRISM_1FB_1EA_0V: hps = &refprism_1fb_1ea_0v; break; case HP_PRISM_1FA_1E_0V: hps = &refprism_1fa_1e_0v; break; case HP_PRISM_2FA_1E_0V: hps = &refprism_2fa_1e_0v; break; case HP_PRISM_1FA_1FB_0E_0V: hps = &refprism_1fa_1fb_0e_0v; break; case HP_PRISM_2FA_1FB_0E_0V: hps = &refprism_2fa_1fb_0e_0v; break; case HP_PRISM_1FA_1FB_1EA_0V: hps = &refprism_1fa_1fb_1ea_0v; break; case HP_PRISM_1FA_1FB_1EB_0V: hps = &refprism_1fa_1fb_1eb_0v; break; case HP_PRISM_2FA_1FB_1EA_0V: hps = &refprism_2fa_1fb_1ea_0v; break; case HP_PRISM_1FB_1EC_0V: hps = &refprism_1fb_1ec_0v; break; case HP_PRISM_1FA_1FB_1EC_0V: hps = &refprism_1fa_1fb_1ec_0v; break; case HP_PRISM_2FA_1FB_1EC_0V: hps = &refprism_2fa_1fb_1ec_0v; break; case HP_PRISM_1FB_2EA_0V: hps = &refprism_1fb_2ea_0v; break; case HP_PRISM_1FA_1FB_2EA_0V: hps = &refprism_1fa_1fb_2ea_0v; break; case HP_PRISM_2FA_1FB_2EA_0V: hps = &refprism_2fa_1fb_2ea_0v; break; case HP_PRISM_1FB_2EB_0V: hps = &refprism_1fb_2eb_0v; break; case HP_PRISM_1FA_1FB_2EB_0V: hps = &refprism_1fa_1fb_2eb_0v; break; case HP_PRISM_1FA_1FB_2EC_0V: hps = &refprism_1fa_1fb_2ec_0v; break; case HP_PRISM_2FA_1FB_2EB_0V: hps = &refprism_2fa_1fb_2eb_0v; break; case HP_PRISM_1FB_3E_0V: hps = &refprism_1fb_3e_0v; break; case HP_PRISM_1FA_1FB_3E_0V: hps = &refprism_1fa_1fb_3e_0v; break; case HP_PRISM_2FA_1FB_3E_0V: hps = &refprism_2fa_1fb_3e_0v; break; case HP_PRISM_2FB_0E_0V: hps = &refprism_2fb_0e_0v; break; case HP_PRISM_1FA_2FB_0E_0V: hps = &refprism_1fa_2fb_0e_0v; break; case HP_PRISM_2FA_2FB_0E_0V: hps = &refprism_2fa_2fb_0e_0v; break; case HP_PRISM_2FB_1EC_0V: hps = &refprism_2fb_1ec_0v; break; case HP_PRISM_1FA_2FB_1EC_0V: hps = &refprism_1fa_2fb_1ec_0v; break; case HP_PRISM_2FA_2FB_1EC_0V: hps = &refprism_2fa_2fb_1ec_0v; break; case HP_PRISM_1FA_2FB_1EB_0V: hps = &refprism_1fa_2fb_1eb_0v; break; case HP_PRISM_2FB_3E_0V: hps = &refprism_2fb_3e_0v; break; case HP_PRISM_1FA_2FB_3E_0V: hps = &refprism_1fa_2fb_3e_0v; break; case HP_PRISM_2FA_2FB_3E_0V: hps = &refprism_2fa_2fb_3e_0v; break; case HP_PRISM_1FA_2E_0V: hps = &refprism_1fa_2e_0v; break; case HP_PRISM_2FA_2E_0V: hps = &refprism_2fa_2e_0v; break; case HP_PRISM_3E_0V: hps = &refprism_3e_0v; break; case HP_PRISM_1FA_3E_0V: hps = &refprism_1fa_3e_0v; break; case HP_PRISM_2FA_3E_0V: hps = &refprism_2fa_3e_0v; break; case HP_PRISM_3FB_0V: hps = &refprism_3fb_0v; break; case HP_PRISM_1FA_3FB_0V: hps = &refprism_1fa_3fb_0v; break; case HP_PRISM_2FA_3FB_0V: hps = &refprism_2fa_3fb_0v; break; // case HP_PRISM_3E_4EH: // hps = &refprism_3e_4eh; break; /*case HP_PRISM_1FB_1EB_0V: hps = &refprism_1fb_1eb_0v; break; case HP_PRISM_2F_0E_0V: hps = &refprism_2f_0e_0v; break; */ case HP_PYRAMID: hps = &refpyramid; break; case HP_PYRAMID_0E_1V: hps = &refpyramid_0e_1v; break; case HP_PYRAMID_EDGES: hps = &refpyramid_edges; break; case HP_PYRAMID_1FB_0E_0V: hps = &refpyramid_1fb_0e_0v; break; case HP_PYRAMID_1FB_0E_1VA: hps = &refpyramid_1fb_0e_1va; break; case HP_HEX: hps = &refhex; break; case HP_HEX_0E_1V: hps = &refhex_0e_1v; break; case HP_HEX_1E_1V: hps = &refhex_1e_1v; break; case HP_HEX_1E_0V: hps = &refhex_1e_0v; break; case HP_HEX_3E_0V: hps = &refhex_3e_0v; break; case HP_HEX_1F_0E_0V: hps = &refhex_1f_0e_0v; break; case HP_HEX_1FA_1FB_0E_0V: hps = &refhex_1fa_1fb_0e_0v; break; case HP_HEX7: hps = &refhex7; break; case HP_HEX7_1FA: hps = &refhex7_1fa; break; case HP_HEX7_1FB: hps = &refhex7_1fb; break; default: { if (GetHPRegistry().count(type)) hps = GetHPRegistry()[type]; else hps = NULL; } } /* if (type != HP_TET_1E_4V && type != HP_TET_1E_2VD) { if (hps->geom == HP_TET) hps = &reftet; if (hps->geom == HP_TRIG) hps = &reftrig; } */ if (!hps) { cout << "Attention hps : hp-refinement not implemented for case " << type << endl; PrintSysError ("hp-refinement not implemented for case ", type); } return hps; } bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & edges, INDEX_2_HASHTABLE & edgepoiclt_dom, NgBitArray & cornerpoint, NgBitArray & edgepoint, INDEX_3_HASHTABLE & faces, INDEX_2_HASHTABLE & face_edges, INDEX_2_HASHTABLE & surf_edges, NgArray & facepoint, int & levels, int & act_ref); bool ClassifyHPElements (Mesh & mesh, NgArray & elements, SplittingType split, int & act_ref, int & levels); void InitHPElements(Mesh & mesh, NgArray & elements) { for(ElementIndex i = 0; i < mesh.GetNE(); i++) { HPRefElement hpel(mesh[i]); hpel.coarse_elnr = i; switch (mesh[i].GetType()) { case PRISM: hpel.type = HP_PRISM; break; case HEX: hpel.type = HP_HEX; break; case TET: hpel.type = HP_TET; break; case PYRAMID: hpel.type = HP_PYRAMID; break; default: cerr << "HPRefElement: illegal elementtype (1) " << mesh[i].GetType() << endl; throw NgException ("HPRefElement: illegal elementtype (1)"); } elements.Append(hpel); } for(SurfaceElementIndex i = 0; i < mesh.GetNSE(); i++) { HPRefElement hpel(mesh[i]); hpel.coarse_elnr = i; switch(mesh[i].GetType()) { case TRIG: hpel.type = HP_TRIG; break; case QUAD: hpel.type = HP_QUAD; break; default: cerr << "HPRefElement: illegal elementtype (1b) " << mesh[i].GetType() << endl; throw NgException ("HPRefElement: illegal elementtype (1b)"); } elements.Append(hpel); } for(SegmentIndex i = 0; i < mesh.GetNSeg(); i++) { Segment & seg = mesh[i]; HPRefElement hpel(mesh[i]); hpel.coarse_elnr = i; hpel.type = HP_SEGM; // hpel.index = seg.edgenr + 10000*seg.si; hpel.index = seg.edgenr; hpel.si = seg.si; /* if(seg.edgenr >= 10000) { throw NgException("assumption that seg.edgenr < 10000 is wrong"); } */ elements.Append(hpel); } } /* ******************************* DoRefinement *************************************** */ void DoRefinement (Mesh & mesh, NgArray & elements, Refinement * ref, double fac1) { elements.SetAllocSize (5 * elements.Size()); INDEX_2_HASHTABLE newpts(elements.Size()+1); INDEX_3_HASHTABLE newfacepts(elements.Size()+1); double fac2 = max(0.001,min(1.0/3,fac1)); // factor for face points PrintMessage(3, " in HP-REFINEMENT with fac1 ", fac1); *testout << " in HP-REFINEMENT with fac1 " << fac1 << endl; int oldelsize = elements.Size(); for (int i = 0; i < oldelsize; i++) { HPRefElement & el = elements[i]; HPRef_Struct * hprs = Get_HPRef_Struct (el.type); if (!hprs) { cout << "Refinementstruct not defined for element " << el.type << endl; continue; } int j = 0; while (hprs->splitedges[j][0]) { INDEX_2 i2(el.pnums[hprs->splitedges[j][0]-1], el.pnums[hprs->splitedges[j][1]-1]); if (fac1 == 0.5) i2.Sort(); if (!newpts.Used (i2)) { Point<3> np; for( int l=0;l<3;l++) np(l) = (1-fac1)*mesh.Point(i2.I1())(l) + fac1 * mesh.Point(i2.I2())(l); int npi = mesh.AddPoint (np); newpts.Set (i2, npi); } j++; } j = 0; if (hprs->splitfaces) while (hprs->splitfaces[j][0]) { INDEX_3 i3(el.pnums[hprs->splitfaces[j][0]-1], el.pnums[hprs->splitfaces[j][1]-1], el.pnums[hprs->splitfaces[j][2]-1]); if (i3.I2() > i3.I3()) Swap (i3.I2(), i3.I3()); if (!newfacepts.Used (i3)) { Point<3> np; for( int l=0;l<3;l++) np(l) = (1-2*fac2)*mesh.Point(i3.I1())(l) + fac2*mesh.Point(i3.I2())(l) + fac2*mesh.Point(i3.I3())(l); int npi = mesh.AddPoint (np); newfacepts.Set (i3, npi); } j++; } } for (int i = 0; i < oldelsize; i++) { HPRefElement el = elements[i]; HPRef_Struct * hprs = Get_HPRef_Struct (el.type); int newlevel = el.levelx; int oldnp = 0; switch (hprs->geom) { case HP_SEGM: oldnp = 2; break; case HP_TRIG: oldnp = 3; break; case HP_QUAD: oldnp = 4; break; case HP_TET: oldnp = 4; break; case HP_PYRAMID: oldnp = 5; break; case HP_PRISM: oldnp = 6; break; case HP_HEX: oldnp = 8; break; case HP_HEX7: oldnp = 7; break; default: cerr << "HPRefElement: illegal type (3) " << hprs->geom << endl; throw NgException ("HPRefElement::SetType: illegal type (3)"); } if (el.type == HP_SEGM || el.type == HP_TRIG || el.type == HP_QUAD || el.type == HP_TET || el.type == HP_PRISM || el.type == HP_HEX || el.type == HP_HEX7 || el.type == HP_PYRAMID) newlevel = el.levelx; if (!hprs) continue; int newpnums[64]; double newparam[64][3]; int j; for (j = 0; j < oldnp; j++) { newpnums[j] = el.pnums[j]; for (int l = 0; l < 3; l++) newparam[j][l] = el.param[j][l]; } // split edges, incl. transferring curvature j = 0; while (hprs->splitedges[j][0]) { INDEX_2 i2(el.pnums[hprs->splitedges[j][0]-1], el.pnums[hprs->splitedges[j][1]-1]); if (fac1 == 0.5) i2.Sort(); int npi = newpts.Get(i2); newpnums[hprs->splitedges[j][2]-1] = npi; for (int l = 0; l < 3; l++) newparam[hprs->splitedges[j][2]-1][l] = (1-fac1) * el.param[hprs->splitedges[j][0]-1][l] + fac1 * el.param[hprs->splitedges[j][1]-1][l]; j++; } // split faces j = 0; if (hprs->splitfaces) while (hprs->splitfaces[j][0]) { INDEX_3 i3(el.pnums[hprs->splitfaces[j][0]-1], el.pnums[hprs->splitfaces[j][1]-1], el.pnums[hprs->splitfaces[j][2]-1]); if (i3.I2() > i3.I3()) Swap (i3.I2(), i3.I3()); int npi = newfacepts.Get(i3); newpnums[hprs->splitfaces[j][3]-1] = npi; for (int l = 0; l < 3; l++) newparam[hprs->splitfaces[j][3]-1][l] = (1-2*fac2) * el.param[hprs->splitfaces[j][0]-1][l] + fac2 * el.param[hprs->splitfaces[j][1]-1][l] + fac2 * el.param[hprs->splitfaces[j][2]-1][l]; j++; } // split elements j = 0; if (hprs->splitelements) while (hprs->splitelements[j][0]) { //int pi1 = el.pnums[hprs->splitelements[j][0]-1]; Point<3> np; for( int l=0;l<3;l++) np(l) = (1-3*fac1)* mesh.Point(el.pnums[hprs->splitelements[j][0]-1])(l) + fac1* mesh.Point(el.pnums[hprs->splitelements[j][1]-1])(l) + fac1* mesh.Point(el.pnums[hprs->splitelements[j][2]-1])(l) + fac1* mesh.Point(el.pnums[hprs->splitelements[j][3]-1])(l); int npi = mesh.AddPoint (np); newpnums[hprs->splitelements[j][4]-1] = npi; for (int l = 0; l < 3; l++) newparam[hprs->splitelements[j][4]-1][l] = (1-3*fac1) * el.param[hprs->splitelements[j][0]-1][l] + fac1 * el.param[hprs->splitelements[j][1]-1][l] + fac1 * el.param[hprs->splitelements[j][2]-1][l] + fac1 * el.param[hprs->splitelements[j][3]-1][l]; j++; } j = 0; /* *testout << " newpnums = "; for (int hi = 0; hi < 64; hi++) *testout << newpnums[hi] << " "; *testout << endl; */ while (hprs->neweltypes[j]) { HPRef_Struct * hprsnew = Get_HPRef_Struct (hprs->neweltypes[j]); HPRefElement newel(el); newel.type = hprs->neweltypes[j]; // newel.index = elements[i].index; // newel.coarse_elnr = elements[i].coarse_elnr; if (newel.type == HP_SEGM || newel.type == HP_TRIG || newel.type == HP_QUAD || newel.type == HP_TET || newel.type == HP_PRISM || newel.type == HP_HEX || newel.type == HP_HEX7 || newel.type == HP_PYRAMID) newel.levelx = newel.levely = newel.levelz = newlevel; else newel.levelx = newel.levely = newel.levelz = newlevel+1; switch(hprsnew->geom) { case HP_SEGM: newel.np=2; break; case HP_QUAD: newel.np=4; break; case HP_TRIG: newel.np=3; break; case HP_HEX: newel.np=8; break; case HP_HEX7: newel.np=7; break; case HP_PRISM: newel.np=6; break; case HP_TET: newel.np=4; break; case HP_PYRAMID: newel.np=5; break; default: throw NgException (string("hprefinement.cpp: illegal type")); } for (int k = 0; k < newel.np; k++) newel.pnums[k] = newpnums[hprs->newels[j][k]-1]; /* *testout << " newel pnums " ; for (int k = 0; k < newel.np; k++) *testout << newel.pnums[k] << "\t"; *testout << endl; */ for (int k = 0; k < newel.np; k++) { for (int l = 0; l < 3; l++) { newel.param[k][l] = newparam[hprs->newels[j][k]-1][l]; // *testout << newel.param[k][l] << " \t "; } // *testout << endl; } if (j == 0) elements[i] = newel; // overwrite old element else elements.Append (newel); j++; } } } /* ************************** DoRefineDummies ******************************** */ void DoRefineDummies (Mesh & mesh, NgArray & elements, Refinement * ref) { int oldelsize = elements.Size(); for (int i = 0; i < oldelsize; i++) { HPRefElement el = elements[i]; HPRef_Struct * hprs = Get_HPRef_Struct (el.type); if (!hprs) continue; if (el.type != HP_DUMMY_QUAD_SINGCORNER && el.type != HP_PYRAMID_EDGES && el.type != HP_PYRAMID_0E_1V && el.type != HP_HEX_0E_1V && el.type != HP_HEX_1E_1V && el.type != HP_HEX_1E_0V && el.type != HP_HEX_3E_0V ) continue; int newlevel = el.levelx; int newpnums[8]; int j; for (j = 0; j < 8; j++) newpnums[j] = el.pnums[j]; double newparam[8][3]; for (j = 0; j < 8; j++) for (int k = 0; k < 3; k++) newparam[j][k] = el.param[j][k]; j = 0; while (hprs->neweltypes[j]) { HPRef_Struct * hprsnew = Get_HPRef_Struct (hprs->neweltypes[j]); HPRefElement newel(el); switch(hprsnew->geom) { case HP_SEGM: newel.np=2; break; case HP_QUAD: newel.np=4; break; case HP_TRIG: newel.np=3; break; case HP_HEX: newel.np=8; break; case HP_PRISM: newel.np=6; break; case HP_TET: newel.np=4; break; case HP_PYRAMID: newel.np=5; break; default: cerr << "HPRefElement: illegal type (4) " << hprsnew->geom << endl; throw NgException ("HPRefElement: illegal type (4)"); } newel.type = hprs->neweltypes[j]; for (int k = 0; k < 8; k++) newel.pnums[k] = newpnums[hprs->newels[j][k]-1]; newel.index = el.index; newel.si = el.si; newel.coarse_elnr = el.coarse_elnr; newel.levelx = newel.levely = newel.levelz = newlevel; for (int k = 0; k < 8; k++) for (int l = 0; l < 3; l++) newel.param[k][l] = newparam[hprs->newels[j][k]-1][l]; if (j == 0) elements[i] = newel; else elements.Append (newel); j++; } } } void SubdivideDegeneratedHexes (Mesh & mesh, NgArray & elements, double fac1) { int oldne = elements.Size(); for (int i = 0; i < oldne; i++) if (Get_HPRef_Struct (elements[i].type)->geom == HP_HEX) { bool common = 0; for (int j = 0; j < 8; j++) for (int k = 0; k < j; k++) if (elements[i].pnums[j] == elements[i].pnums[k]) common = 1; if (common) { cout << " Degenerate Hex found " << endl; *testout << " Degenerate Hex found " << endl; HPRefElement el = elements[i]; HPRefElement newel = el; Point<3> center(0,0,0); double newparam[3] = { 0, 0, 0 }; for (int j = 0; j < 8; j++) { center += 0.125 * Vec<3>(mesh[el.pnums[j]]); // 0.125 originates form 8 points not from fac1; for (int l = 0; l < 3; l++) newparam[l] += 0.125 * el.param[j][l]; } int npi = mesh.AddPoint (center); const ELEMENT_FACE * faces = MeshTopology::GetFaces1 (HEX); for (int j = 0; j < 6; j++) { NgArray pts; for (int k = 0; k < 4; k++) { bool same = 0; for (int l = 0; l < pts.Size(); l++) if (el.pnums[pts[l]] == el.pnums[faces[j][k]-1]) same = 1; if (!same) pts.Append (faces[j][k]-1); } if (pts.Size() == 3) // TrigFace -> TET { for (int k = 0; k < 3; k++) { newel.pnums[k] = el.pnums[pts[2-k]]; for (int l = 0; l < 3; l++) newel.param[k][l] = el.param[pts[2-k]][l]; } newel.pnums[3] = npi; for (int l = 0; l < 3; l++) newel.param[3][l] = newparam[l]; newel.type = HP_TET; newel.np = 4; } else { for (int k = 0; k < 4; k++) { newel.pnums[k] = el.pnums[pts[3-k]]; for (int l = 0; l < 3; l++) newel.param[k][l] = el.param[pts[3-k]][l]; } newel.pnums[4] = npi; for (int l = 0; l < 3; l++) newel.param[4][l] = newparam[l]; newel.type = HP_PYRAMID; newel.np = 5; } if (j == 0) elements[i] = newel; else elements.Append (newel); } /* const ELEMENT_EDGE * edges = MeshTopology::GetEdges (HEX); for(int k=0;k<12;k++) { int e[2]; for(int l=0;l<2;l++) e[l] = edges[k][l]-1; if(el.PNum(e[0]+1)!=el.PNum(e[1]+1)) { newel.SetType(HP_SEGM); for(int l=0;l<2;l++) { newel.pnums[0] = el.PNum(e[l]+1); newel.pnums[1] = npi; for(int j=0;j<3;j++) { // newel.param[0][j] = el.param[e[l]][j]; // newel.param[1][j] = newparam[j]; } elements.Append(newel); } newel.SetType(HP_TRIG); newel.pnums[0] = el.PNum(e[0]+1); newel.pnums[1] = el.PNum(e[1]+1); newel.pnums[2] = npi; *testout << "DEGHEX TRIG :: newpnums " << newel.pnums[0] << "\t" << newel.pnums[1] << "\t" << newel.pnums[2] << endl; cout << "DEGHEX TRIG :: newpnums " << newel.pnums[0] << "\t" << newel.pnums[1] << "\t" << newel.pnums[2] << endl; for(int j=0;j<3;j++) { // newel.param[0][j] = el.param[e[0]][j]; // newel.param[1][j] = el.param[e[1]][j]; // newel.param[2][j] = newparam[j]; } elements.Append(newel); } }*/ } } } void CalcStatistics (NgArray & elements) { return; #ifdef ABC int i, p; int nsegm = 0, ntrig = 0, nquad = 0; int nhex = 0, nprism = 0, npyramid = 0, ntet = 0; int maxlevel = 0; for (i = 1; i <= elements.Size(); i++) { const HPRefElement & el = elements.Get(i); maxlevel = max2 (el.level, maxlevel); switch (Get_HPRef_Struct (el.type)->geom) { case HP_SEGM: { nsegm++; break; } case HP_TRIG: { ntrig ++; break; } case HP_QUAD: { nquad++; break; } case HP_TET: { ntet++; break; } case HP_PRISM: { nprism++; break; } case HP_PYRAMID: { npyramid++; break; } case HP_HEX: { nhex++; break; } default: { cerr << "statistics error, unknown element type" << endl; } } } cout << "level = " << maxlevel << endl; cout << "nsegm = " << nsegm << endl; cout << "ntrig = " << ntrig << ", nquad = " << nquad << endl; cout << "ntet = " << ntet << ", npyr = " << npyramid << ", nprism = " << nprism << ", nhex = " << nhex << endl; return; double memcost = 0, cpucost = 0; for (p = 1; p <= 20; p++) { memcost = (ntet + nprism + nhex) * pow (static_cast(p), 6.0); cpucost = (ntet + nprism + nhex) * pow (static_cast(p), 9.0); cout << "costs for p = " << p << ": mem = " << memcost << ", cpu = " << cpucost << endl; } double memcosttet = 0; double memcostprism = 0; double memcosthex = 0; double memcostsctet = 0; double memcostscprism = 0; double memcostschex = 0; double cpucosttet = 0; double cpucostprism = 0; double cpucosthex = 0; for (i = 1; i <= elements.Size(); i++) { const HPRefElement & el = elements.Get(i); switch (el.type) { case HP_TET: case HP_TET_0E_1V: case HP_TET_1E_0V: case HP_TET_1E_1VA: { int p1 = maxlevel - el.level + 1; (*testout) << "p1 = " << p1 << ", P1^6 = " << pow (static_cast(p1), 6.0) << " (p1-3)^6 = " << pow ( static_cast(max2(p1-3, 0)), 6.0) << " p1^3 = " << pow ( static_cast(p1), 3.0) << " (p1-3)^3 = " << pow ( static_cast(p1-3), 3.0) << " [p1^3-(p1-3)^3]^2 = " << sqr (pow (static_cast(p1),3.0) - pow ( static_cast(p1-3), 3.0)) << endl; p1 /= 2 +1; memcosttet += pow (static_cast(p1), 6.0); memcostsctet += pow (static_cast(p1), 6.0) - pow ( static_cast(max2(p1-3, 1)), 6.0); cpucosttet += pow (static_cast(p1), 9.0); break; } case HP_PRISM: case HP_PRISM_SINGEDGE: { int p1 = maxlevel - el.level + 1; p1 /= 2 +1; memcostprism += pow (static_cast(p1), 6.0); memcostscprism += pow (static_cast(p1), 6.0) - pow ( static_cast(max2(p1-3, 1)), 6.0); cpucostprism += pow (static_cast(p1), 9.0); break; } case HP_HEX: { int p1 = maxlevel - el.level + 1; int p2 = maxlevel; p1 /= 2 +1; p2 /= 2 +1; memcosthex += pow (static_cast(p1), 4.0) * pow (static_cast(p2), 2.0); memcostschex += pow (static_cast(p1), 6.0) - pow ( static_cast(max2(p1-2, 0)), 6.0); cpucosthex += pow (static_cast(p1), 6.0) * pow (static_cast(p2), 3.0); break; } default: ; } } cout << "TET: hp-memcost = " << memcosttet << ", scmemcost = " << memcostsctet << ", cpucost = " << cpucosttet << endl; cout << "PRI: hp-memcost = " << memcostprism << ", scmemcost = " << memcostscprism << ", cpucost = " << cpucostprism << endl; cout << "HEX: hp-memcost = " << memcosthex << ", scmemcost = " << memcostschex << ", cpucost = " << cpucosthex << endl; #endif } void ReorderPoints (Mesh & mesh, NgArray & hpelements) { NgArray map (mesh.GetNP()); for (int i = 1; i <= mesh.GetNP(); i++) map[i] = i; int nwrong(0), nright(0); for (int k = 0; k < 5; k++) { nwrong = nright = 0; for (int i = 0; i < hpelements.Size(); i++) { const HPRefElement & hpel = hpelements[i]; if (Get_HPRef_Struct (hpel.type) -> geom == HP_PRISM) { int minbot = 0, mintop = 0; for (int j = 0; j < 3; j++) { if (map[hpel.pnums[j]] < map[hpel.pnums[minbot]]) minbot = j; if (map[hpel.pnums[j+3]] < map[hpel.pnums[mintop+3]]) mintop = j; } if (minbot != mintop) nwrong++; else nright++; if (minbot != mintop) { if (map[hpel.pnums[minbot]] < map[hpel.pnums[mintop+3]]) swap (map[hpel.pnums[3+minbot]], map[hpel.pnums[3+mintop]]); else swap (map[hpel.pnums[minbot]], map[hpel.pnums[mintop]]); } } } // cout << nwrong << " wrong prisms, " << nright << " right prisms" << endl; } PrintMessage(3, nwrong, " wrong prisms, ", nright, " right prisms"); NgArray hpts(mesh.GetNP()); for (int i = 1; i <= mesh.GetNP(); i++) hpts[map[i]] = mesh.Point(i); for (int i = 1; i <= mesh.GetNP(); i++) mesh.Point(i) = hpts[i]; for (int i = 0; i < hpelements.Size(); i++) { HPRefElement & hpel = hpelements[i]; for (int j = 0; j < hpel.np; j++) hpel.pnums[j] = map[hpel.pnums[j]]; } } /* ***************************** HPRefinement ********************************** */ void HPRefinement (Mesh & mesh, Refinement * ref, SplittingType split, int levels, double fac1, bool setorders, bool reflevels) { PrintMessage (1, "HP Refinement called, levels = ", levels); // NgLock mem_lock (mem_mutex,1); mesh.coarsemesh = make_unique(); *mesh.coarsemesh = mesh; // #ifdef CURVEDELEMS_NEW const_cast (mesh.coarsemesh->GetCurvedElements() ). BuildCurvedElements (ref, mesh.GetCurvedElements().GetOrder()); // #endif // delete mesh.hpelements; mesh.hpelements = make_unique>(); NgArray & hpelements = *mesh.hpelements; InitHPElements(mesh,hpelements); NgArray nplevel; nplevel.Append (mesh.GetNP()); int act_ref=1; bool sing = ClassifyHPElements (mesh,hpelements, split, act_ref, levels); sing = true; // iterate at least once while(sing) { PrintMessage(3, " Start new hp-refinement: step ", act_ref); DoRefinement (mesh, hpelements, ref, fac1); DoRefineDummies (mesh, hpelements, ref); nplevel.Append (mesh.GetNP()); CalcStatistics (hpelements); SubdivideDegeneratedHexes (mesh, hpelements,fac1); ReorderPoints (mesh, hpelements); mesh.ClearSegments(); mesh.ClearSurfaceElements(); mesh.ClearVolumeElements(); for (int i = 0; i < hpelements.Size(); i++) { HPRefElement & hpel = hpelements[i]; if (Get_HPRef_Struct (hpel.type)) switch (Get_HPRef_Struct (hpel.type) -> geom) { case HP_SEGM: { Segment seg; seg[0] = hpel.pnums[0]; seg[1] = hpel.pnums[1]; // NOTE: only for less than 10000 elements (HACK) !!! // seg.edgenr = hpel.index % 10000; // seg.si = hpel.index / 10000; seg.edgenr = hpel.index; seg.si = hpel.si; /* seg.epgeominfo[0].dist = hpel.param[0][0]; // he: war hpel.param[0][0] seg.epgeominfo[1].dist = hpel.param[1][0]; // he: war hpel.param[1][0] */ const Segment & coarseseg = mesh.coarsemesh->LineSegment(hpel.coarse_elnr+1); double d1 = coarseseg.epgeominfo[0].dist; double d2 = coarseseg.epgeominfo[1].dist; // seg.epgeominfo[0].dist = hpel.param[0][0]; // he: war hpel.param[0][0] // seg.epgeominfo[1].dist = hpel.param[1][0]; // he: war hpel.param[1][0] seg.epgeominfo[0].dist = d1 + hpel.param[0][0] * (d2-d1); // JS, June 08 seg.epgeominfo[1].dist = d1 + hpel.param[1][0] * (d2-d1); seg.epgeominfo[0].edgenr = seg.edgenr; seg.epgeominfo[1].edgenr = seg.edgenr; seg.domin = hpel.domin; seg.domout=hpel.domout; // he: needed for segments! seg.hp_elnr = i; seg.singedge_left = hpel.singedge_left; seg.singedge_right = hpel.singedge_right; seg.SetCurved (coarseseg.IsCurved()); mesh.AddSegment (seg); break; } case HP_TRIG: case HP_QUAD: { Element2d el(hpel.np); for(int j=0;j geom)); } } PrintMessage(5, " Start with Update Topology "); mesh.UpdateTopology(); PrintMessage(5, " Mesh Update Topology done "); act_ref++; sing = ClassifyHPElements(mesh,hpelements, split, act_ref, levels); } PrintMessage(3, " HP-Refinement done with ", --act_ref, " refinement steps."); if(act_ref>=1) { for(ElementIndex i=0;i v(hpel.param[edges[j][0]-1][0]-hpel.param[edges[j][1]-1][0], hpel.param[edges[j][0]-1][1]-hpel.param[edges[j][1]-1][1], hpel.param[edges[j][0]-1][2]-hpel.param[edges[j][1]-1][2]); dist[edge_dir[j]] = max(v.Length(),dist[edge_dir[j]]); } int refi[3]; for(int j=0;j<3;j++) refi[j] = int(max(double(floor(log(dist[ord_dir[j]]/sqrt(2.))/log(fac1))),0.)); // cout << " ref " << refi[0] << "\t" << refi[1] << "\t" << refi[2] << endl; // cout << " order " << act_ref +1 - refi[0] << "\t" << act_ref +1 - refi[1] << "\t" << act_ref +1 - refi[2] << endl; if(setorders) mesh[i].SetOrder(act_ref+1-refi[0],act_ref+1-refi[1],act_ref+1-refi[2]); } for(SurfaceElementIndex i=0;i v(hpel.param[edges[j][0]-1][0]-hpel.param[edges[j][1]-1][0], hpel.param[edges[j][0]-1][1]-hpel.param[edges[j][1]-1][1], hpel.param[edges[j][0]-1][2]-hpel.param[edges[j][1]-1][2]); dist[edge_dir[j]] = max(v.Length(),dist[edge_dir[j]]); } int refi[3]; for(int j=0;j<3;j++) refi[j] = int(max(double(floor(log(dist[ord_dir[j]]/sqrt(2.))/log(fac1))),0.)); if(setorders) mesh[i].SetOrder(act_ref+1-refi[0],act_ref+1-refi[1],act_ref+1-refi[2]); // cout << " ref " << refi[0] << "\t" << refi[1] << endl; // cout << " order " << act_ref +1 - refi[0] << "\t" << act_ref +1 - refi[1] << endl; } } } bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & edges, INDEX_2_HASHTABLE & edgepoint_dom, NgBitArray & cornerpoint, NgBitArray & edgepoint, INDEX_3_HASHTABLE & faces, INDEX_2_HASHTABLE & face_edges, INDEX_2_HASHTABLE & surf_edges, NgArray & facepoint, int & levels, int & act_ref) { bool sing = 0; if (mesh.GetDimension() == 3) { /* // check, if point has as least 3 different surfs: NgArray surfonpoint(mesh.GetNP()); surfonpoint = INDEX_3(0,0,0); for (SurfaceElementIndex sei = 0; sei < mesh.GetNSE(); sei++) { const Element2d & el = mesh[sei]; int ind = el.GetIndex(); for (int j = 0; j < el.GetNP(); j++) { INDEX_3 & i3 = surfonpoint[el[j]]; if (ind != i3.I1() && ind != i3.I2() && ind != i3.I3()) { i3.I1() = i3.I2(); i3.I2() = i3.I3(); i3.I3() = ind; } } } for (int i = 1; i <= mesh.GetNP(); i++) if (surfonpoint.Get(i).I1()) cornerpoint.Set(i); */ cornerpoint.Clear(); for (int i = 1; i <= mesh.GetNP(); i++) { if (mesh.Point(i).Singularity() * levels >= act_ref) { cornerpoint.Set(i); sing = 1; } } // cout << endl; for (int i = 1; i <= mesh.GetNSeg(); i++) if (mesh.LineSegment(i).singedge_left * levels >= act_ref) { INDEX_2 i2 (mesh.LineSegment(i)[0], mesh.LineSegment(i)[1]); /* // before edges.Set (i2, 1); i2.Sort(); INDEX_2 i2s(i2.I2(), i2.I1()); edges.Set (i2s, 1); */ edges.Set (i2, 1); INDEX_2 i2s(i2.I2(), i2.I1()); edges.Set (i2s, 1); edgepoint.Set (i2.I1()); edgepoint.Set (i2.I2()); sing = 1; } // if 2 adjacent edges of an element are singular, the // common point must be a singular point for (int i = 1; i <= mesh.GetNE(); i++) { const Element & el = mesh.VolumeElement(i); const ELEMENT_EDGE * eledges = MeshTopology::GetEdges1 (el.GetType()); int nedges = MeshTopology::GetNEdges (el.GetType()); for (int j = 0; j < nedges; j++) for (int k = 0; k < nedges; k++) if (j != k) { INDEX_2 ej(el.PNum(eledges[j][0]), el.PNum(eledges[j][1])); ej.Sort(); INDEX_2 ek(el.PNum(eledges[k][0]), el.PNum(eledges[k][1])); ek.Sort(); if (edges.Used(ej) && edges.Used(ek)) { if (ej.I1() == ek.I1()) cornerpoint.Set (ek.I1()); if (ej.I1() == ek.I2()) cornerpoint.Set (ek.I2()); if (ej.I2() == ek.I1()) cornerpoint.Set (ek.I1()); if (ej.I2() == ek.I2()) cornerpoint.Set (ek.I2()); } } } edgepoint.Or (cornerpoint); (*testout) << "cornerpoint = " << endl << cornerpoint << endl; (*testout) << "edgepoint = " << endl << edgepoint << endl; facepoint = 0; for (SurfaceElementIndex sei = 0; sei < mesh.GetNSE(); sei++) { const Element2d & el = mesh[sei]; const FaceDescriptor & fd = mesh.GetFaceDescriptor (el.GetIndex()); int domnr = 0; if (fd.DomainInSingular() * levels < act_ref && fd.DomainOutSingular() * levels < act_ref) { domnr=0; continue;} if (fd.DomainInSingular() * levels >= act_ref) { domnr = fd.DomainIn(); sing = 1; } if (fd.DomainOutSingular() * levels >= act_ref) { domnr = fd.DomainOut(); sing = 1; } if (fd.DomainInSingular() * levels >= act_ref && fd.DomainOutSingular() * levels >= act_ref) { domnr = -1; sing = 1; } INDEX_3 i3; if (el.GetNP() == 3) i3 = INDEX_3::Sort (el[0], el[1], el[2]); else { INDEX_4 i4 (el[0], el[1], el[2], el[3]); i4.Sort(); i3 = INDEX_3(i4.I1(), i4.I2(), i4.I3()); } faces.Set (i3, domnr); *testout << "set face " << i3 << ", domnr = " << domnr << endl; for (int j = 0; j < el.GetNP(); j++) { face_edges.Set (INDEX_2::Sort (el[j], el[(j+1)%el.GetNP()]), domnr); surf_edges.Set (INDEX_2::Sort (el[j], el[(j+1)%el.GetNP()]), fd.SurfNr()+1); facepoint[el[j]] = domnr; } } // (*testout) << "singular edges = " << edges << endl; // (*testout) << "singular faces = " << faces << endl; // (*testout) << "singular faces_edges = " << face_edges << endl; } else { // 2D case // check, if point has as least 3 different surfs: NgArray surfonpoint(mesh.GetNP()); for (int i = 1; i <= mesh.GetNP(); i++) surfonpoint.Elem(i) = INDEX_3(0,0,0); for (int i = 1; i <= mesh.GetNSeg(); i++) { const Segment & seg = mesh.LineSegment(i); int ind = seg.edgenr; if (seg.singedge_left * levels >= act_ref) { INDEX_2 i2 = INDEX_2::Sort(mesh.LineSegment(i)[0], mesh.LineSegment(i)[1]); edges.Set(i2,1); edgepoint.Set(i2.I1()); edgepoint.Set(i2.I2()); *testout << " singleft " << endl; *testout << " mesh.LineSegment(i).domout " << mesh.LineSegment(i).domout << endl; *testout << " mesh.LineSegment(i).domin " << mesh.LineSegment(i).domin << endl; edgepoint_dom.Set (INDEX_2(mesh.LineSegment(i).domin, i2.I1()), 1); edgepoint_dom.Set (INDEX_2(mesh.LineSegment(i).domin, i2.I2()), 1); sing = 1; } if (seg.singedge_right * levels >= act_ref) { INDEX_2 i2 = INDEX_2::Sort(mesh.LineSegment(i)[1], mesh.LineSegment(i)[0]); edges.Set (i2, 1); edgepoint.Set(i2.I1()); edgepoint.Set(i2.I2()); *testout << " singright " << endl; *testout << " mesh.LineSegment(i).domout " << mesh.LineSegment(i).domout << endl; *testout << " mesh.LineSegment(i).domin " << mesh.LineSegment(i).domin << endl; edgepoint_dom.Set (INDEX_2(mesh.LineSegment(i).domout, i2.I1()), 1); edgepoint_dom.Set (INDEX_2(mesh.LineSegment(i).domout, i2.I2()), 1); sing = 1; } // (*testout) << "seg = " << ind << ", " << seg[0] << "-" << seg[1] << endl; if (seg.singedge_left * levels >= act_ref || seg.singedge_right* levels >= act_ref) { for (int j = 0; j < 2; j++) { int pi = (j == 0) ? seg[0] : seg[1]; INDEX_3 & i3 = surfonpoint.Elem(pi); if (ind != i3.I1() && ind != i3.I2()) { i3.I1() = i3.I2(); i3.I2() = ind; } } } } for (int i = 1; i <= mesh.GetNP(); i++) { // mark points for refinement that are in corners between two anisotropic edges if (surfonpoint.Get(i).I1()) { // cornerpoint.Set(i); // disabled by JS, Aug 2009 edgepoint.Set(i); } // mark points for refinement that are explicitly specified in input file if (mesh.Point(i).Singularity()*levels >= act_ref) { cornerpoint.Set(i); edgepoint.Set(i); sing = 1; } } edgepoint.Or (cornerpoint); (*testout) << "2d sing edges: " << endl << edges << endl; (*testout) << "2d cornerpoints: " << endl << cornerpoint << endl << "2d edgepoints: " << endl << edgepoint << endl; facepoint = 0; } if (!sing) PrintMessage(3, "PrepareElements no more to do for actual refinement ", act_ref); return(sing); } bool ClassifyHPElements (Mesh & mesh, NgArray & elements, SplittingType split, int & act_ref, int & levels) { INDEX_2_HASHTABLE edges(mesh.GetNSeg()+1); NgBitArray edgepoint(mesh.GetNP()); INDEX_2_HASHTABLE edgepoint_dom(mesh.GetNSeg()+1); edgepoint.Clear(); NgBitArray cornerpoint(mesh.GetNP()); cornerpoint.Clear(); // value = nr > 0 ... refine elements in domain nr // value = -1 ..... refine elements in any domain INDEX_3_HASHTABLE faces(mesh.GetNSE()+1); INDEX_2_HASHTABLE face_edges(mesh.GetNSE()+1); INDEX_2_HASHTABLE surf_edges(mesh.GetNSE()+1); NgArray facepoint(mesh.GetNP()); bool sing = CheckSingularities(mesh, edges, edgepoint_dom, cornerpoint, edgepoint, faces, face_edges, surf_edges, facepoint, levels, act_ref); if (act_ref == 1 && split == SPLIT_ALFELD) sing = true; if (act_ref == 1 && split == SPLIT_POWELL) sing = true; if(sing==0) return(sing); int cnt_undef = 0, cnt_nonimplement = 0; NgArray misses(10000); misses = 0; (*testout) << "edgepoint_dom = " << endl << edgepoint_dom << endl; for( int i = 0; igeom) { case HP_TET: { hpel.type = ClassifyTet(hpel, edges, edgepoint_dom, cornerpoint, edgepoint, faces,face_edges, surf_edges, facepoint); /* // if (i != 182) if ( (!hpel.PNums().Contains(40)) || (!hpel.PNums().Contains(41)) ) hpel.type = HP_NONETET; else cout << "type = " << hpel.type << endl; */ break; } case HP_PRISM: { hpel.type = ClassifyPrism(hpel, edges, edgepoint_dom, cornerpoint, edgepoint, faces, face_edges, surf_edges, facepoint); break; } case HP_HEX7: { hpel.type = ClassifyHex7(hpel, edges, edgepoint_dom, cornerpoint, edgepoint, faces, face_edges, surf_edges, facepoint); break; } case HP_HEX: { hpel.type = ClassifyHex(hpel, edges, edgepoint_dom, cornerpoint, edgepoint, faces, face_edges, surf_edges, facepoint); break; } case HP_TRIG: { int dim = mesh.GetDimension(); const FaceDescriptor & fd = mesh.GetFaceDescriptor (hpel.GetIndex()); if (split == SPLIT_HP) hpel.type = ClassifyTrig(hpel, edges, edgepoint_dom, cornerpoint, edgepoint, faces, face_edges, surf_edges, facepoint, dim, fd); else if (split == SPLIT_ALFELD) hpel.type = HP_TRIG_ALFELD; else if (split == SPLIT_POWELL) hpel.type = HP_TRIG_POWELL; dd = 2; break; } case HP_QUAD: { int dim = mesh.GetDimension(); const FaceDescriptor & fd = mesh.GetFaceDescriptor (hpel.GetIndex()); hpel.type = ClassifyQuad(hpel, edges, edgepoint_dom, cornerpoint, edgepoint, faces, face_edges, surf_edges, facepoint, dim, fd); dd = 2; break; } case HP_SEGM: { if (split == SPLIT_HP) hpel.type = ClassifySegm(hpel, edges, edgepoint_dom, cornerpoint, edgepoint, faces, face_edges, surf_edges, facepoint); else if (split == SPLIT_POWELL) hpel.type = HP_SEGM_SINGCORNERL; dd = 1; break; } case HP_PYRAMID: { hpel.type = ClassifyPyramid(hpel, edges, edgepoint_dom, cornerpoint, edgepoint, faces, face_edges, surf_edges, facepoint); break; } default: { cout << "illegal element type for hp-prepare elements " << hpel.type << endl; throw NgException ("hprefinement.cpp: don't know how to set parameters"); } } if(hpel.type == HP_NONE) cnt_undef++; //else //cout << "elem " << i << " classified type " << hpel.type << endl; if (!Get_HPRef_Struct (hpel.type)) { (*testout) << "hp-element-type " << hpel.type << " not implemented " << endl; (*testout) << " elType " << hprs->geom << endl; (cout) << " elType " << hprs->geom << endl; cnt_nonimplement++; misses[hpel.type]++; } for(int j=0; j