mirror of
https://github.com/NGSolve/netgen.git
synced 2024-09-23 02:02:49 +05:00
1963 lines
53 KiB
C++
1963 lines
53 KiB
C++
#include <mystdlib.h>
|
|
#include "meshing.hpp"
|
|
#include "hprefinement.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) :
|
|
np(el.GetNV()), index(el.GetIndex()), levelx(0), levely(0), levelz(0), type(HP_NONE), domin(-1), domout(-1) //domin,out for segements
|
|
{
|
|
//Reset();
|
|
for (int i=0; i<np ; i++)
|
|
pnums[i] = el[i];
|
|
|
|
const Point3d * points =
|
|
MeshTopology :: GetVertices (el.GetType());
|
|
for(int i=0;i<np;i++)
|
|
for(int l=0;l<3;l++)
|
|
param[i][l] = points[i].X(l+1);
|
|
}
|
|
|
|
|
|
HPRefElement :: HPRefElement(Element2d & el) :
|
|
levelx(0), levely(0), levelz(0), type(HP_NONE), index(el.GetIndex()), np(el.GetNV()), domin(-1), domout(-1) //domin,out for segements
|
|
{
|
|
//Reset();
|
|
|
|
for (int i=0; i<np ; i++)
|
|
pnums[i] = el[i];
|
|
|
|
const Point3d * points =
|
|
MeshTopology :: GetVertices (el.GetType());
|
|
for(int i=0;i<np;i++)
|
|
for(int l=0;l<3;l++)
|
|
param[i][l] = points[i].X(l+1);
|
|
}
|
|
|
|
HPRefElement :: HPRefElement(Segment & el) :
|
|
levelx(0), levely(0), levelz(0), type(HP_NONE), np(2), domin(el.domin), domout(el.domout), singedge_left(el.singedge_left), singedge_right(el.singedge_right)
|
|
{
|
|
//Reset();
|
|
for (int i=0; i<np ; i++)
|
|
pnums[i] = el[i];
|
|
const Point3d * points =
|
|
MeshTopology :: GetVertices (SEGMENT);
|
|
for(int i=0;i<np;i++)
|
|
for(int l=0;l<3;l++)
|
|
param[i][l] = points[i].X(l+1);
|
|
|
|
/*
|
|
for (int i=0; i<np; i++)
|
|
{
|
|
param[i][0] = i;
|
|
param[i][1] = -1; param[i][2] = -1;
|
|
}
|
|
*/
|
|
}
|
|
|
|
HPRefElement :: HPRefElement(HPRefElement & el) :
|
|
np(el.np), levelx(el.levelx), levely(el.levely), levelz(el.levelz), type(el.type), domin(el.domin), domout(el.domout), index(el.index), coarse_elnr(el.coarse_elnr), singedge_left(el.singedge_left), singedge_right(el.singedge_right)
|
|
|
|
{
|
|
//Reset();
|
|
for (int i=0; i<np ; i++)
|
|
{
|
|
pnums[i] = el[i];
|
|
for(int l=0; l<3; l++) param[i][l] = el.param[i][l];
|
|
}
|
|
}
|
|
|
|
void HPRefElement :: SetType (HPREF_ELEMENT_TYPE t)
|
|
{
|
|
type = t;
|
|
switch(type)
|
|
{
|
|
case HP_SEGM: np=2; break;
|
|
case HP_TRIG: np=3; break;
|
|
case HP_QUAD: np=4; break;
|
|
case HP_TET: np=4; break;
|
|
case HP_PRISM: np=6; break;
|
|
case HP_PYRAMID: np=5; break;
|
|
case HP_HEX: np=8; break;
|
|
|
|
default:
|
|
cerr << "HPRefElement: illegal type " << type << endl;
|
|
throw NgException ("HPRefElement::SetType: illegal type");
|
|
}
|
|
|
|
for(int k = 0; k < 8;k++)
|
|
{
|
|
pnums[k]=0;
|
|
for(int l = 0; l < 3; l++)
|
|
param[k][l]=0.;
|
|
}
|
|
}
|
|
|
|
|
|
HPRef_Struct * Get_HPRef_Struct (HPREF_ELEMENT_TYPE type)
|
|
{
|
|
HPRef_Struct * hps = NULL;
|
|
|
|
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_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_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_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_1EA_0V:
|
|
hps = &reftet_1f_1ea_0v; break;
|
|
case HP_TET_1F_1EB_0V:
|
|
hps = &reftet_1f_1eb_0v; break;
|
|
case HP_TET_2F_0E_0V:
|
|
hps = &reftet_2f_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_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;
|
|
|
|
|
|
|
|
default:
|
|
{
|
|
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<int> & edges, INDEX_2_HASHTABLE<int> & edgepoiclt_dom,
|
|
NgBitArray & cornerpoint, NgBitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges,
|
|
INDEX_2_HASHTABLE<int> & surf_edges, NgArray<int, PointIndex::BASE> & facepoint, int & levels, int & act_ref);
|
|
|
|
bool ClassifyHPElements (Mesh & mesh, NgArray<HPRefElement> & elements, int & act_ref, int & levels);
|
|
|
|
|
|
void InitHPElements(Mesh & mesh, NgArray<HPRefElement> & 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;
|
|
if(seg.edgenr >= 10000)
|
|
{
|
|
throw NgException("assumption that seg.edgenr < 10000 is wrong");
|
|
}
|
|
elements.Append(hpel);
|
|
}
|
|
}
|
|
|
|
|
|
|
|
/* ******************************* DoRefinement *************************************** */
|
|
void DoRefinement (Mesh & mesh, NgArray<HPRefElement> & elements,
|
|
Refinement * ref, double fac1)
|
|
{
|
|
elements.SetAllocSize (5 * elements.Size());
|
|
INDEX_2_HASHTABLE<int> newpts(elements.Size()+1);
|
|
INDEX_3_HASHTABLE<int> newfacepts(elements.Size()+1);
|
|
|
|
// prepare new points
|
|
|
|
fac1 = max(0.001,min(0.33,fac1));
|
|
cout << " in HP-REFINEMENT with fac1 " << fac1 << endl;
|
|
*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 (!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*fac1)*mesh.Point(i3.I1())(l)
|
|
+ fac1*mesh.Point(i3.I2())(l) + fac1*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;
|
|
|
|
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_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]);
|
|
|
|
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*fac1) * el.param[hprs->splitfaces[j][0]-1][l] +
|
|
fac1 * el.param[hprs->splitfaces[j][1]-1][l] +
|
|
fac1 * 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_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_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<HPRefElement> & 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.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<HPRefElement> & 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<int> 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<HPRefElement> & 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<double>(p), 6.0);
|
|
cpucost = (ntet + nprism + nhex) * pow (static_cast<double>(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<double>(p1), 6.0)
|
|
<< " (p1-3)^6 = " << pow ( static_cast<double>(max2(p1-3, 0)), 6.0)
|
|
<< " p1^3 = " << pow ( static_cast<double>(p1), 3.0)
|
|
<< " (p1-3)^3 = " << pow ( static_cast<double>(p1-3), 3.0)
|
|
<< " [p1^3-(p1-3)^3]^2 = " << sqr (pow (static_cast<double>(p1),3.0) - pow ( static_cast<double>(p1-3), 3.0))
|
|
<< endl;
|
|
|
|
p1 /= 2 +1;
|
|
memcosttet += pow (static_cast<double>(p1), 6.0);
|
|
memcostsctet += pow (static_cast<double>(p1), 6.0) - pow ( static_cast<double>(max2(p1-3, 1)), 6.0);
|
|
cpucosttet += pow (static_cast<double>(p1), 9.0);
|
|
break;
|
|
}
|
|
case HP_PRISM:
|
|
case HP_PRISM_SINGEDGE:
|
|
{
|
|
int p1 = maxlevel - el.level + 1;
|
|
p1 /= 2 +1;
|
|
memcostprism += pow (static_cast<double>(p1), 6.0);
|
|
memcostscprism += pow (static_cast<double>(p1), 6.0) - pow ( static_cast<double>(max2(p1-3, 1)), 6.0);
|
|
cpucostprism += pow (static_cast<double>(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<double>(p1), 4.0) * pow (static_cast<double>(p2), 2.0);
|
|
memcostschex += pow (static_cast<double>(p1), 6.0) - pow ( static_cast<double>(max2(p1-2, 0)), 6.0);
|
|
cpucosthex += pow (static_cast<double>(p1), 6.0) * pow (static_cast<double>(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<HPRefElement> & hpelements)
|
|
{
|
|
NgArray<int, 1> 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;
|
|
}
|
|
|
|
cout << nwrong << " wrong prisms, " << nright << " right prisms" << endl;
|
|
|
|
|
|
NgArray<MeshPoint, 1> 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, 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>();
|
|
*mesh.coarsemesh = mesh;
|
|
|
|
// #ifdef CURVEDELEMS_NEW
|
|
const_cast<CurvedElements&> (mesh.coarsemesh->GetCurvedElements() ).
|
|
BuildCurvedElements (ref, mesh.GetCurvedElements().GetOrder());
|
|
// #endif
|
|
|
|
|
|
// delete mesh.hpelements;
|
|
mesh.hpelements = make_unique<NgArray<HPRefElement>>();
|
|
|
|
NgArray<HPRefElement> & hpelements = *mesh.hpelements;
|
|
|
|
InitHPElements(mesh,hpelements);
|
|
|
|
NgArray<int> nplevel;
|
|
nplevel.Append (mesh.GetNP());
|
|
|
|
int act_ref=1;
|
|
bool sing = ClassifyHPElements (mesh,hpelements, act_ref, levels);
|
|
|
|
sing = true; // iterate at least once
|
|
while(sing)
|
|
{
|
|
cout << " Start new hp-refinement: step " << act_ref << endl;
|
|
|
|
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.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;
|
|
mesh.AddSegment (seg);
|
|
break;
|
|
}
|
|
|
|
case HP_TRIG:
|
|
case HP_QUAD:
|
|
{
|
|
Element2d el(hpel.np);
|
|
for(int j=0;j<hpel.np;j++)
|
|
el.PNum(j+1) = hpel.pnums[j];
|
|
el.hp_elnr = i;
|
|
el.SetIndex(hpel.index);
|
|
if(setorders)
|
|
el.SetOrder(act_ref+1,act_ref+1,0);
|
|
mesh.AddSurfaceElement(el);
|
|
break;
|
|
}
|
|
case HP_HEX:
|
|
case HP_TET:
|
|
case HP_PRISM:
|
|
case HP_PYRAMID:
|
|
{
|
|
Element el(hpel.np);
|
|
for(int j=0;j<hpel.np;j++)
|
|
el.PNum(j+1) = hpel.pnums[j];
|
|
el.SetIndex(hpel.index);
|
|
el.hp_elnr = i;
|
|
if(setorders)
|
|
el.SetOrder(act_ref+1,act_ref+1,act_ref+1);
|
|
mesh.AddVolumeElement(el);
|
|
break;
|
|
}
|
|
|
|
default:
|
|
PrintSysError ("hpref, backconversion failed for element ",
|
|
int(Get_HPRef_Struct (hpel.type) -> geom));
|
|
}
|
|
}
|
|
cout << " Start with Update Topology " << endl;
|
|
mesh.UpdateTopology();
|
|
cout << " Mesh Update Topology done " << endl;
|
|
|
|
act_ref++;
|
|
|
|
sing = ClassifyHPElements(mesh,hpelements, act_ref, levels);
|
|
}
|
|
|
|
cout << " HP-Refinement done with " << --act_ref << " refinement steps." << endl;
|
|
|
|
if(act_ref>=1)
|
|
{
|
|
for(ElementIndex i=0;i<mesh.GetNE(); i++)
|
|
{
|
|
// Element el = mesh[i] ;
|
|
HPRefElement & hpel = hpelements[mesh[i].hp_elnr];
|
|
const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (mesh[i].GetType());
|
|
double dist[3] = {0,0,0};
|
|
int ord_dir[3] = {0,0,0};
|
|
int edge_dir[12] = {0,0,0,0,0,0,0,0,0,0,0,0};
|
|
int ned = 4;
|
|
|
|
switch (mesh[i].GetType())
|
|
{
|
|
case TET:
|
|
/* cout << " TET " ;
|
|
for(int k=0;k<4;k++) cout << el[k] << "\t" ;
|
|
cout << endl; */
|
|
break;
|
|
case PRISM:
|
|
/* cout << " PRISM " ;
|
|
for(int k=0;k<6;k++) cout << el[k] << "\t" ;
|
|
cout << endl; */
|
|
for(int l=6;l<9;l++) edge_dir[l] = 2;
|
|
ord_dir[2] = 2;
|
|
ned = 9;
|
|
break;
|
|
case HEX:
|
|
/* cout << " HEX " ;
|
|
for(int k=0;k<8;k++) cout << el[k] << "\t" ;
|
|
cout << endl; */
|
|
for(int l=8;l<12; l++) edge_dir[l] = 2;
|
|
edge_dir[2] = edge_dir[3] = edge_dir[6] = edge_dir[7] = 1;
|
|
ord_dir[1] = 1;
|
|
ord_dir[2] = 2;
|
|
ned = 12;
|
|
break;
|
|
case PYRAMID:
|
|
/* cout << " PYRAMID " ;
|
|
for(int k=0;k<5;k++) cout << el[k] << "\t" ;
|
|
cout << endl; */
|
|
for(int l=4;l<8;l++) edge_dir[l] = 2;
|
|
edge_dir[2] = edge_dir[3] = 1;
|
|
ord_dir[1] = 1;
|
|
ord_dir[2] = 2;
|
|
ned = 8;
|
|
break;
|
|
|
|
|
|
default:
|
|
cerr << "HPRefElement: illegal elementtype (2) " << mesh[i].GetType() << endl;
|
|
throw NgException ("HPRefElement: illegal elementtype (2)");
|
|
|
|
}
|
|
|
|
for (int j=0;j<ned;j++)
|
|
{
|
|
|
|
Vec<3> 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<mesh.GetNSE(); i++)
|
|
{
|
|
// Element2d el = mesh[i] ;
|
|
HPRefElement & hpel = hpelements[mesh[i].hp_elnr];
|
|
const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (mesh[i].GetType());
|
|
double dist[3] = {0,0,0};
|
|
int ord_dir[3] = {0,0,0};
|
|
int edge_dir[4] = {0,0,0,0} ;
|
|
int ned = 3;
|
|
|
|
if(mesh[i].GetType() == QUAD)
|
|
{
|
|
/* cout << " QUAD " ;
|
|
for(int k=0;k<4;k++) cout << el[k] << "\t" ;
|
|
cout << endl; */
|
|
|
|
edge_dir[2] = edge_dir[3] = 1;
|
|
ord_dir[1] = 1;
|
|
ned = 4;
|
|
}
|
|
/* else
|
|
{
|
|
cout << " TRIG " ;
|
|
for(int k=0;k<3;k++) cout << el[k] << "\t" ;
|
|
cout << endl;
|
|
} */
|
|
|
|
for (int j=0;j<ned;j++)
|
|
{
|
|
Vec<3> 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<int> & edges, INDEX_2_HASHTABLE<int> & edgepoint_dom,
|
|
NgBitArray & cornerpoint, NgBitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges,
|
|
INDEX_2_HASHTABLE<int> & surf_edges, NgArray<int, PointIndex::BASE> & facepoint, int & levels, int & act_ref)
|
|
{
|
|
bool sing = 0;
|
|
if (mesh.GetDimension() == 3)
|
|
{
|
|
/*
|
|
// check, if point has as least 3 different surfs:
|
|
|
|
NgArray<INDEX_3, PointIndex::BASE> 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
|
|
// commen 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);
|
|
|
|
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 faces = " << faces << endl;
|
|
(*testout) << "singular faces_edges = " << face_edges << endl;
|
|
}
|
|
else
|
|
{
|
|
// 2D case
|
|
|
|
// check, if point has as least 3 different surfs:
|
|
NgArray<INDEX_3, PointIndex::BASE> 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 (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 (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)
|
|
cout << "PrepareElements no more to do for actual refinement " << act_ref << endl;
|
|
|
|
return(sing);
|
|
}
|
|
|
|
|
|
|
|
bool ClassifyHPElements (Mesh & mesh, NgArray<HPRefElement> & elements, int & act_ref, int & levels)
|
|
{
|
|
INDEX_2_HASHTABLE<int> edges(mesh.GetNSeg()+1);
|
|
NgBitArray edgepoint(mesh.GetNP());
|
|
INDEX_2_HASHTABLE<int> 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<int> faces(mesh.GetNSE()+1);
|
|
INDEX_2_HASHTABLE<int> face_edges(mesh.GetNSE()+1);
|
|
INDEX_2_HASHTABLE<int> surf_edges(mesh.GetNSE()+1);
|
|
NgArray<int, PointIndex::BASE> facepoint(mesh.GetNP());
|
|
|
|
bool sing = CheckSingularities(mesh, edges, edgepoint_dom,
|
|
cornerpoint, edgepoint, faces, face_edges,
|
|
surf_edges, facepoint, levels, act_ref);
|
|
|
|
if(sing==0) return(sing);
|
|
|
|
int cnt_undef = 0, cnt_nonimplement = 0;
|
|
NgArray<int> misses(10000);
|
|
misses = 0;
|
|
|
|
(*testout) << "edgepoint_dom = " << endl << edgepoint_dom << endl;
|
|
|
|
for( int i = 0; i<elements.Size(); i++)
|
|
{
|
|
// *testout << "classify element " << i << endl;
|
|
|
|
HPRefElement & hpel = elements[i];
|
|
HPRef_Struct * hprs = Get_HPRef_Struct (hpel.type);
|
|
HPRefElement old_el = elements[i];
|
|
int dd=3;
|
|
|
|
|
|
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))
|
|
continue;
|
|
|
|
sing = 1;
|
|
switch (hprs->geom)
|
|
{
|
|
case HP_TET:
|
|
{
|
|
hpel.type = ClassifyTet(hpel, edges, edgepoint_dom, cornerpoint, edgepoint, faces,face_edges, surf_edges, facepoint);
|
|
break;
|
|
}
|
|
case HP_PRISM:
|
|
{
|
|
hpel.type = ClassifyPrism(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());
|
|
|
|
hpel.type = ClassifyTrig(hpel, edges, edgepoint_dom, cornerpoint, edgepoint,
|
|
faces, face_edges, surf_edges, facepoint, dim, fd);
|
|
|
|
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:
|
|
{
|
|
hpel.type = ClassifySegm(hpel, edges, edgepoint_dom, cornerpoint, edgepoint,
|
|
faces, face_edges, surf_edges, facepoint);
|
|
dd = 1;
|
|
break;
|
|
}
|
|
case HP_PYRAMID:
|
|
{
|
|
hpel.type = ClassifyPyramid(hpel, edges, edgepoint_dom, cornerpoint, edgepoint, faces,
|
|
face_edges, surf_edges, facepoint);
|
|
|
|
cout << " ** Pyramid classified " << hpel.type << endl;
|
|
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<hpel.np; j++)
|
|
{
|
|
for( int k=0; k<hpel.np; k++)
|
|
if(hpel[j] == old_el.pnums[k])
|
|
{
|
|
for(int l=0;l<dd;l++)
|
|
hpel.param[j][l] = old_el.param[k][l];
|
|
break;
|
|
}
|
|
}
|
|
|
|
}
|
|
|
|
|
|
cout << "undefined elements update classification: " << cnt_undef << endl;
|
|
cout << "non-implemented in update classification: " << cnt_nonimplement << endl;
|
|
|
|
for (int i = 0; i < misses.Size(); i++)
|
|
if (misses[i])
|
|
cout << " in update classification missing case " << i << " occurred " << misses[i] << " times" << endl;
|
|
|
|
return(sing);
|
|
}
|
|
}
|
|
|