ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2024-12-26 13:50:33 +05:00
Code Issues Packages Projects Releases Wiki Activity

edge/face topology

Browse Source
This commit is contained in:
Joachim Schoeberl 2014-01-07 10:42:30 +00:00
parent 081e9f129e
commit 3d8e1d407f
1 changed files with 159 additions and 102 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

171
libsrc/visualization/vssolution.cpp
View File

@ -2393,7 +2393,7 @@ namespace netgen
MyMPI_Bcast (comp); MyMPI_Bcast (comp);
#endif #endif
double val; // double val;
// bool considerElem; // bool considerElem;
bool hasit = false; bool hasit = false;
@ -2413,32 +2413,28 @@ namespace netgen
NgProfiler::RegionTimer reg1 (timer1); NgProfiler::RegionTimer reg1 (timer1);
int ne = mesh->GetNE(); int ne = mesh->GetNE();
double hmax = maxv;
double hmin = minv;
#pragma omp parallel for reduction (max : hmax) reduction (min : hmin) double hminv = numeric_limits<double>::max();
double hmaxv = -numeric_limits<double>::max();
bool hhasit = false;
#pragma omp parallel for reduction (max : hmaxv) reduction (min : hminv) reduction (|| : hhasit)
for (int i = 0; i < ne; i++) for (int i = 0; i < ne; i++)
{ {
double val;
bool considerElem = GetValue (sol, i, 0.333, 0.333, 0.333, comp, val); bool considerElem = GetValue (sol, i, 0.333, 0.333, 0.333, comp, val);
if (considerElem) if (considerElem)
{ {
/* if (val > hmaxv) hmaxv = val;
#pragma omp critical(getminmaxvol) if (val < hminv) hminv = val;
{ hhasit = true;
if (val > maxv || !hasit) maxv = val;
if (val < minv || !hasit) minv = val;
hasit = true;
}
*/
if (val > hmax) hmax = val;
if (val < hmin) hmin = val;
} }
} }
maxv = hmax; minv = min(minv, hminv);
minv = hmin; maxv = max(maxv, hmaxv);
// cout << "maxv = " << maxv << " =?= " << hmax << endl; hasit |= hhasit;
} }
if (sol->draw_surface) if (sol->draw_surface)
{ {
NgProfiler::RegionTimer reg2 (timer2); NgProfiler::RegionTimer reg2 (timer2);
@ -2447,6 +2443,7 @@ namespace netgen
for (int i = 0; i < nse; i++) for (int i = 0; i < nse; i++)
{ {
ELEMENT_TYPE type = mesh->SurfaceElement(i+1).GetType(); ELEMENT_TYPE type = mesh->SurfaceElement(i+1).GetType();
double val;
bool considerElem = (type == QUAD) bool considerElem = (type == QUAD)
? GetSurfValue (sol, i, -1, 0.5, 0.5, comp, val) ? GetSurfValue (sol, i, -1, 0.5, 0.5, comp, val)
: GetSurfValue (sol, i, -1, 0.3333333, 0.3333333, comp, val); : GetSurfValue (sol, i, -1, 0.3333333, 0.3333333, comp, val);
@ -2461,6 +2458,7 @@ namespace netgen
} }
} }
if (minv == maxv) maxv = minv+1e-6; if (minv == maxv) maxv = minv+1e-6;
if (!hasit) { minv = 0; maxv = 1; }
#ifdef PARALLEL #ifdef PARALLEL
if ((ntasks > 1) && (id == 0)) if ((ntasks > 1) && (id == 0))
@ -3656,14 +3654,18 @@ namespace netgen
void VisualSceneSolution :: GetClippingPlaneTrigs (Array<ClipPlaneTrig> & trigs, void VisualSceneSolution :: GetClippingPlaneTrigs (Array<ClipPlaneTrig> & trigs,
Array<ClipPlanePoint> & pts) Array<ClipPlanePoint> & pts)
{ {
static int timer_vals = NgProfiler::CreateTimer ("ClipPlaneTrigs - vertex values");
static int timer1 = NgProfiler::CreateTimer ("ClipPlaneTrigs1"); static int timer1 = NgProfiler::CreateTimer ("ClipPlaneTrigs1");
static int timer2 = NgProfiler::CreateTimer ("ClipPlaneTrigs2"); static int timer1a = NgProfiler::CreateTimer ("ClipPlaneTrigs1a");
// static int timer2 = NgProfiler::CreateTimer ("ClipPlaneTrigs2");
static int timer3 = NgProfiler::CreateTimer ("ClipPlaneTrigs3"); static int timer3 = NgProfiler::CreateTimer ("ClipPlaneTrigs3");
static int timer4 = NgProfiler::CreateTimer ("ClipPlaneTrigs4"); static int timer4 = NgProfiler::CreateTimer ("ClipPlaneTrigs4");
static int timer4b = NgProfiler::CreateTimer ("ClipPlaneTrigs4b");
NgProfiler::RegionTimer reg1 (timer1); NgProfiler::RegionTimer reg1 (timer1);
int ne = mesh->GetNE(); int ne = mesh->GetNE();
const int edgei[6][2] = const int edgei[6][2] =
@ -3688,11 +3690,30 @@ namespace netgen
Array<Point<3> > locgrid(n3); Array<Point<3> > locgrid(n3);
Array<Mat<3,3> > trans(n3); Array<Mat<3,3> > trans(n3);
Array<double> val(n3); Array<double> val(n3);
Array<bool> locposval(n3);
Array<int> compress(n3); Array<int> compress(n3);
NgProfiler::StartTimer (timer_vals);
Array<double,PointIndex::BASE> vertval(mesh->GetNV());
Array<bool,PointIndex::BASE> posval(mesh->GetNV());
for (PointIndex pi = vertval.Begin(); pi < vertval.End(); pi++)
{
Point<3> vert = (*mesh)[pi];
vertval[pi] =
vert(0) * clipplane[0] +
vert(1) * clipplane[1] +
vert(2) * clipplane[2] +
clipplane[3];
posval[pi] = vertval[pi] > 0;
}
NgProfiler::StopTimer (timer_vals);
INDEX_2_CLOSED_HASHTABLE<int> edges(8*n3); // point nr of edge
for (ElementIndex ei = 0; ei < ne; ei++) for (ElementIndex ei = 0; ei < ne; ei++)
{ {
// NgProfiler::RegionTimer reg1a (timer1a);
int first_point_of_element = pts.Size(); int first_point_of_element = pts.Size();
locgrid.SetSize(n3); locgrid.SetSize(n3);
@ -3707,7 +3728,26 @@ namespace netgen
int ii = 0; int ii = 0;
int cnt_valid = 0; int cnt_valid = 0;
NgProfiler::StartTimer (timer2); // NgProfiler::StartTimer (timer2);
if (!mesh->GetCurvedElements().IsElementHighOrder(ei))
{
bool has_pos = 0, has_neg = 0;
for (int i = 0; i < el.GetNP(); i++)
if (posval[el[i]] > 0)
has_pos = 1;
else
has_neg = 1;
if (!has_pos || !has_neg)
{
// NgProfiler::StopTimer (timer2);
continue;
}
}
if (type == TET || type == TET10) if (type == TET || type == TET10)
@ -3770,16 +3810,18 @@ namespace netgen
locgrid.SetSize(cnt_valid); locgrid.SetSize(cnt_valid);
NgProfiler::StopTimer (timer2); // NgProfiler::StopTimer (timer2);
NgProfiler::RegionTimer reg4(timer4); // NgProfiler::RegionTimer reg4(timer4);
if (mesh->GetCurvedElements().IsHighOrder()) if (mesh->GetCurvedElements().IsHighOrder())
{ {
NgProfiler::RegionTimer reg4(timer4);
mesh->GetCurvedElements(). mesh->GetCurvedElements().
CalcMultiPointElementTransformation (&locgrid, ei, &grid, 0); CalcMultiPointElementTransformation (&locgrid, ei, &grid, 0);
} }
else else
{ {
NgProfiler::RegionTimer reg4(timer4b);
Vector shape(el.GetNP()); Vector shape(el.GetNP());
MatrixFixWidth<3> pointmat(el.GetNP()); MatrixFixWidth<3> pointmat(el.GetNP());
@ -3803,7 +3845,7 @@ namespace netgen
NgProfiler::RegionTimer reg3(timer3); NgProfiler::RegionTimer reg3(timer3);
bool has_pos = 0, has_neg = 0; bool has_pos = false, all_pos = true;
for (int i = 0; i < cnt_valid; i++) for (int i = 0; i < cnt_valid; i++)
{ {
@ -3813,14 +3855,17 @@ namespace netgen
grid[i](2) * clipplane[2] + grid[i](2) * clipplane[2] +
clipplane[3]; clipplane[3];
if (val[i] > 0) locposval[i] = val[i] > 0;
has_pos = 1; has_pos |= locposval[i];
else all_pos &= locposval[i];
has_neg = 1;
// if (val[i] > 0) has_pos = 1; else has_neg = 1;
} }
if (!has_pos || !has_neg) continue; // if (!has_pos || !has_neg) continue;
if (!has_pos || all_pos) continue;
edges.DeleteData();
for (int ix = 0; ix < n; ix++) for (int ix = 0; ix < n; ix++)
for (int iy = 0; iy < n; iy++) for (int iy = 0; iy < n; iy++)
@ -3834,6 +3879,16 @@ namespace netgen
for (int j = 0; j < 8; j++) for (int j = 0; j < 8; j++)
pi[j] = compress[pi[j]]; pi[j] = compress[pi[j]];
bool has_pos = false, all_pos = true;
for (int j = 0; j < 8; j++)
if (pi[j] != -1)
{
has_pos |= locposval[pi[j]];
all_pos &= locposval[pi[j]];
}
if (!has_pos || all_pos) continue;
const int tets[6][4] = const int tets[6][4] =
{ { 1, 2, 4, 5 }, { { 1, 2, 4, 5 },
{ 4, 5, 2, 8 }, { 4, 5, 2, 8 },
@ -3848,11 +3903,19 @@ namespace netgen
for (int k = 0; k < 4; k++) for (int k = 0; k < 4; k++)
teti[k] = pi[tets[ii][k]-1]; teti[k] = pi[tets[ii][k]-1];
bool is_valid = 1; bool is_valid = true;
for (int j = 0; j < 4; j++) for (int j = 0; j < 4; j++)
if (teti[j] == -1) is_valid = 0; is_valid &= (teti[j] != -1);
if (!is_valid) continue; if (!is_valid) continue;
bool has_pos = false, all_pos = true;
for (int j = 0; j < 4; j++)
{
has_pos |= locposval[teti[j]];
all_pos &= locposval[teti[j]];
}
if (!has_pos || all_pos) continue;
for (int j = 0; j < 4; j++) for (int j = 0; j < 4; j++)
nodevali[j] = val[teti[j]]; nodevali[j] = val[teti[j]];
@ -3861,15 +3924,8 @@ namespace netgen
{ {
int lpi1 = edgei[j][0]; int lpi1 = edgei[j][0];
int lpi2 = edgei[j][1]; int lpi2 = edgei[j][1];
if ( (nodevali[lpi1] > 0) != if ( (nodevali[lpi1] > 0) != (nodevali[lpi2] > 0) )
(nodevali[lpi2] > 0) )
{ {
edgelam[j] = nodevali[lpi2] / (nodevali[lpi2] - nodevali[lpi1]);
Point<3> p1 = grid[teti[lpi1]];
Point<3> p2 = grid[teti[lpi2]];
edgep[j] = p1 + (1-edgelam[j]) * (p2-p1);
cntce++; cntce++;
cpe3 = cpe2; cpe3 = cpe2;
cpe2 = cpe1; cpe2 = cpe1;
@ -3888,33 +3944,34 @@ namespace netgen
case 1: ednr = cpe2; break; case 1: ednr = cpe2; break;
case 2: ednr = cpe3; break; case 2: ednr = cpe3; break;
} }
// cpt.points[k].p = edgep[ednr];
int pi1 = edgei[ednr][0]; int pi1 = edgei[ednr][0];
int pi2 = edgei[ednr][1]; int pi2 = edgei[ednr][1];
Point<3> p1 = locgrid[teti[pi1]];
Point<3> p2 = locgrid[teti[pi2]];
// cpt.points[k].lami = p2 + edgelam[ednr] * (p1-p2);
ClipPlanePoint cppt;
cppt.elnr = ei;
cppt.p = edgep[ednr];
cppt.lami = p2 + edgelam[ednr] * (p1-p2);
int pnr = -1; int pnr = -1;
for (int l = first_point_of_element; l < pts.Size(); l++) INDEX_2 pair (teti[pi1], teti[pi2]);
if (fabs (cppt.lami(0)-pts[l].lami(0)) < 1e-8 && pair.Sort();
fabs (cppt.lami(1)-pts[l].lami(1)) < 1e-8 && if (edges.Used(pair))
fabs (cppt.lami(2)-pts[l].lami(2)) < 1e-8) pnr = edges.Get(pair);
else
{ {
pnr = l; ClipPlanePoint cppt;
break; cppt.elnr = ei;
}
edgelam[ednr] = nodevali[pi2] / (nodevali[pi2] - nodevali[pi1]);
Point<3> gp1 = grid[teti[pi1]];
Point<3> gp2 = grid[teti[pi2]];
cppt.p = gp2 + edgelam[ednr] * (gp1-gp2);
Point<3> p1 = locgrid[teti[pi1]];
Point<3> p2 = locgrid[teti[pi2]];
cppt.lami = p2 + edgelam[ednr] * (p1-p2);
if (pnr == -1)
pnr = pts.Append (cppt)-1; pnr = pts.Append (cppt)-1;
edges.Set (pair, pnr);
}
cpt.points[k].pnr = pnr; cpt.points[k].pnr = pnr;
cpt.points[k].locpnr = pnr-first_point_of_element; cpt.points[k].locpnr = pnr-first_point_of_element;
@ -3930,7 +3987,7 @@ namespace netgen
else else
{ // other elements not supported (JS, June 2007) { // other elements not supported (JS, June 2007)
return; continue; // return;
} }
} }