Alefeld split hack

This commit is contained in:
Joachim Schoeberl 2023-07-12 10:50:03 -07:00
parent ce8ec2eb7f
commit 82e88f3afb
3 changed files with 67 additions and 16 deletions

View File

@ -774,3 +774,47 @@ HPRef_Struct reftrig_3singedges =
reftrig_3singedges_newelstypes, reftrig_3singedges_newelstypes,
reftrig_3singedges_newels reftrig_3singedges_newels
}; };
// HP_TRIG_3SINGEDGES
int reftrig_Alefeld_splitedges[][3] =
{
{ 0, 0, 0 }
};
int reftrig_Alefeld_splitfaces[][4] =
{
{ 1, 2, 3, 4 },
{ 0, 0, 0, 0 }
};
HPREF_ELEMENT_TYPE reftrig_Alefeld_newelstypes[] =
{
HP_TRIG, HP_TRIG, HP_TRIG,
HP_NONE,
};
int reftrig_Alefeld_newels[][8] =
{
{ 1, 2, 4 },
{ 2, 3, 4 },
{ 3, 1, 4 },
};
HPRef_Struct reftrig_Alefeld =
{
HP_TRIG,
reftrig_Alefeld_splitedges,
reftrig_Alefeld_splitfaces,
0,
reftrig_Alefeld_newelstypes,
reftrig_Alefeld_newels
};

View File

@ -174,7 +174,10 @@ namespace netgen
case HP_TRIG_3SINGEDGES: case HP_TRIG_3SINGEDGES:
hps = &reftrig_3singedges; break; hps = &reftrig_3singedges; break;
case HP_TRIG_ALEFELD:
hps = &reftrig_Alefeld; break;
case HP_QUAD: case HP_QUAD:
hps = &refquad; break; hps = &refquad; break;
case HP_DUMMY_QUAD_SINGCORNER: case HP_DUMMY_QUAD_SINGCORNER:
@ -1338,7 +1341,7 @@ namespace netgen
sing = true; // iterate at least once sing = true; // iterate at least once
while(sing) while(sing)
{ {
PrintMessage(3, " Start new hp-refinement: step ", act_ref); PrintMessage(0, " Start new hp-refinement: step ", act_ref);
DoRefinement (mesh, hpelements, ref, fac1); DoRefinement (mesh, hpelements, ref, fac1);
DoRefineDummies (mesh, hpelements, ref); DoRefineDummies (mesh, hpelements, ref);
@ -1823,6 +1826,7 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HAS
bool ClassifyHPElements (Mesh & mesh, NgArray<HPRefElement> & elements, int & act_ref, int & levels) bool ClassifyHPElements (Mesh & mesh, NgArray<HPRefElement> & elements, int & act_ref, int & levels)
{ {
cout << "in classify" << endl;
INDEX_2_HASHTABLE<int> edges(mesh.GetNSeg()+1); INDEX_2_HASHTABLE<int> edges(mesh.GetNSeg()+1);
NgBitArray edgepoint(mesh.GetNP()); NgBitArray edgepoint(mesh.GetNP());
INDEX_2_HASHTABLE<int> edgepoint_dom(mesh.GetNSeg()+1); INDEX_2_HASHTABLE<int> edgepoint_dom(mesh.GetNSeg()+1);
@ -1841,7 +1845,9 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HAS
bool sing = CheckSingularities(mesh, edges, edgepoint_dom, bool sing = CheckSingularities(mesh, edges, edgepoint_dom,
cornerpoint, edgepoint, faces, face_edges, cornerpoint, edgepoint, faces, face_edges,
surf_edges, facepoint, levels, act_ref); surf_edges, facepoint, levels, act_ref);
if (act_ref == 1)
sing = true; // Joachim
if(sing==0) return(sing); if(sing==0) return(sing);
int cnt_undef = 0, cnt_nonimplement = 0; int cnt_undef = 0, cnt_nonimplement = 0;
@ -1859,12 +1865,11 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HAS
HPRefElement old_el = elements[i]; HPRefElement old_el = elements[i];
int dd=3; int dd=3;
if(act_ref !=1 && (hpel.type == HP_HEX || hpel.type == HP_PRISM || hpel.type == HP_TET if(act_ref !=1 && (hpel.type == HP_HEX || hpel.type == HP_PRISM || hpel.type == HP_TET
|| hpel.type == HP_PYRAMID || hpel.type == HP_QUAD || hpel.type == HP_TRIG || hpel.type == HP_SEGM)) || hpel.type == HP_PYRAMID || hpel.type == HP_QUAD || hpel.type == HP_TRIG || hpel.type == HP_SEGM))
continue; continue;
sing = 1; sing = 1;
switch (hprs->geom) switch (hprs->geom)
{ {
case HP_TET: case HP_TET:
@ -1887,23 +1892,25 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HAS
break; break;
} }
case HP_TRIG: case HP_TRIG:
{ {
int dim = mesh.GetDimension(); int dim = mesh.GetDimension();
const FaceDescriptor & fd = mesh.GetFaceDescriptor (hpel.GetIndex()); const FaceDescriptor & fd = mesh.GetFaceDescriptor (hpel.GetIndex());
hpel.type = HP_TRIG_ALEFELD;
/*
hpel.type = ClassifyTrig(hpel, edges, edgepoint_dom, cornerpoint, edgepoint, hpel.type = ClassifyTrig(hpel, edges, edgepoint_dom, cornerpoint, edgepoint,
faces, face_edges, surf_edges, facepoint, dim, fd); faces, face_edges, surf_edges, facepoint, dim, fd);
*/
dd = 2; dd = 2;
continue;
break; break;
} }
case HP_QUAD: case HP_QUAD:
{ {
int dim = mesh.GetDimension(); int dim = mesh.GetDimension();
const FaceDescriptor & fd = mesh.GetFaceDescriptor (hpel.GetIndex()); const FaceDescriptor & fd = mesh.GetFaceDescriptor (hpel.GetIndex());
hpel.type = ClassifyQuad(hpel, edges, edgepoint_dom, cornerpoint, edgepoint, hpel.type = ClassifyQuad(hpel, edges, edgepoint_dom, cornerpoint, edgepoint,
faces, face_edges, surf_edges, facepoint, dim, fd); faces, face_edges, surf_edges, facepoint, dim, fd);
dd = 2; dd = 2;
break; break;
} }
@ -1928,6 +1935,8 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HAS
throw NgException ("hprefinement.cpp: don't know how to set parameters"); throw NgException ("hprefinement.cpp: don't know how to set parameters");
} }
} }
continue;
if(hpel.type == HP_NONE) if(hpel.type == HP_NONE)
cnt_undef++; cnt_undef++;
@ -1935,8 +1944,6 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HAS
//else //else
//cout << "elem " << i << " classified type " << hpel.type << endl; //cout << "elem " << i << " classified type " << hpel.type << endl;
if (!Get_HPRef_Struct (hpel.type)) if (!Get_HPRef_Struct (hpel.type))
{ {
(*testout) << "hp-element-type " << hpel.type << " not implemented " << endl; (*testout) << "hp-element-type " << hpel.type << " not implemented " << endl;
@ -1959,8 +1966,7 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE<int> & edges, INDEX_2_HAS
} }
} }
PrintMessage(3, "undefined elements update classification: ", cnt_undef); PrintMessage(3, "undefined elements update classification: ", cnt_undef);
PrintMessage(3, "non-implemented in update classification: ", cnt_nonimplement); PrintMessage(3, "non-implemented in update classification: ", cnt_nonimplement);

View File

@ -41,6 +41,8 @@ enum HPREF_ELEMENT_TYPE {
HP_TRIG_SINGEDGES23, HP_TRIG_SINGEDGES23,
HP_TRIG_3SINGEDGES = 40, HP_TRIG_3SINGEDGES = 40,
HP_TRIG_ALEFELD,
HP_QUAD = 50, HP_QUAD = 50,
HP_QUAD_SINGCORNER, HP_QUAD_SINGCORNER,
HP_DUMMY_QUAD_SINGCORNER, HP_DUMMY_QUAD_SINGCORNER,
@ -311,7 +313,6 @@ public:
}; };
DLL_HEADER extern void HPRefinement (Mesh & mesh, Refinement * ref, int levels, DLL_HEADER extern void HPRefinement (Mesh & mesh, Refinement * ref, int levels,
double fac1=0.125, bool setorders=true, bool ref_level = false); double fac1=0.125, bool setorders=true, bool ref_level = false);