diff --git a/CHANGELOG b/CHANGELOG index 568fdb14..e42bf778 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -14,6 +14,8 @@ Version 4.2.1 (development) - Significantly improved memory usage. +- Add support to visualize solutions on 1D elements embedded in 2D and 3D. + Version 4.2 released on May 23, 2022 ==================================== diff --git a/glvis.cpp b/glvis.cpp index 90c05396..57774f4a 100644 --- a/glvis.cpp +++ b/glvis.cpp @@ -1148,8 +1148,8 @@ int main (int argc, char *argv[]) const char *font_name = string_default; int portnum = 19916; int multisample = GetMultisample(); - double line_width = 1.0; - double ms_line_width = gl3::LINE_WIDTH_AA; + double line_width = GetLineWidth(); + double ms_line_width = GetLineWidthMS(); int geom_ref_type = Quadrature1D::ClosedUniform; bool legacy_gl_ctx = false; bool enable_hidpi = true; diff --git a/lib/openglvis.cpp b/lib/openglvis.cpp index 15b6c2b4..8a626217 100644 --- a/lib/openglvis.cpp +++ b/lib/openglvis.cpp @@ -303,6 +303,25 @@ ::DrawQuad(gl3::GlDrawable& buff, {fpts[3], fnorm, texcoord[3]}); } + +void VisualizationScene +::DrawLine(gl3::GlDrawable& buff, + const double (&pts)[4][3], const double (&cv)[4], + const double minv, const double maxv) +{ + float texcoord[2]; + std::array fpts[2]; + for (int i = 0; i < 2; i++) + { + texcoord[i] = palette.GetColorCoord(cv[i], minv, maxv); + fpts[i] = {(float) pts[i][0], (float) pts[i][1], (float) pts[i][2]}; + } + buff.addLine( + {fpts[0], texcoord[0]}, + {fpts[1], texcoord[1]} + ); +} + void VisualizationScene ::DrawPatch(gl3::GlDrawable& drawable, const DenseMatrix &pts, Vector &vals, DenseMatrix &normals, @@ -376,41 +395,65 @@ ::DrawPatch(gl3::GlDrawable& drawable, const DenseMatrix &pts, Vector &vals, } else { - if (n == 3) + if (n == 2) + { + poly.glBegin(GL_LINES); + } + else if (n == 3) { poly.glBegin(GL_TRIANGLES); } - else + else if (n == 4) { poly.glBegin(GL_QUADS); } + else // We should never hit this case + { + mfem_error("VisualizationScene::DrawPatch() : unknown input geometry."); + } for (int i = 0; i < ind.Size(); i += n) { - int j; - if (n == 3) - j = Compute3DUnitNormal(&pts(0, ind[i]), &pts(0, ind[i+1]), - &pts(0, ind[i+2]), na); - else - j = Compute3DUnitNormal(&pts(0, ind[i]), &pts(0, ind[i+1]), - &pts(0, ind[i+2]), &pts(0, ind[i+3]), na); - if (j == 0) + if (n == 2) // Lines { - if (normals_opt == 0) + // Draw solution line + for (int j=0; j<2; j++) { - poly.glNormal3dv(na); - for ( ; j < n; j++) - { - MySetColor(poly, vals(ind[i+j]), minv, maxv); - poly.glVertex3dv(&pts(0, ind[i+j])); - } + MySetColor(poly, vals(ind[i+j]), minv, maxv); + poly.glVertex3d(pts(0, ind[i+j]), pts(1, ind[i+j]), pts(2, ind[i+j])); + } + } + else // Quads and triangles + { + int j; + if (n == 3) + { + j = Compute3DUnitNormal(&pts(0, ind[i]), &pts(0, ind[i+1]), + &pts(0, ind[i+2]), na); } else { - poly.glNormal3d(-na[0], -na[1], -na[2]); - for (j = n-1; j >= 0; j--) + j = Compute3DUnitNormal(&pts(0, ind[i]), &pts(0, ind[i+1]), + &pts(0, ind[i+2]), &pts(0, ind[i+3]), na); + } + if (j == 0) + { + if (normals_opt == 0) + { + poly.glNormal3dv(na); + for ( ; j < n; j++) + { + MySetColor(poly, vals(ind[i+j]), minv, maxv); + poly.glVertex3dv(&pts(0, ind[i+j])); + } + } + else { - MySetColor(poly, vals(ind[i+j]), minv, maxv); - poly.glVertex3dv(&pts(0, ind[i+j])); + poly.glNormal3d(-na[0], -na[1], -na[2]); + for (j = n-1; j >= 0; j--) + { + MySetColor(poly, vals(ind[i+j]), minv, maxv); + poly.glVertex3dv(&pts(0, ind[i+j])); + } } } } diff --git a/lib/openglvis.hpp b/lib/openglvis.hpp index 676e2f4d..c42d94a9 100644 --- a/lib/openglvis.hpp +++ b/lib/openglvis.hpp @@ -125,6 +125,10 @@ class VisualizationScene const double (&pts)[4][3], const double (&cv)[4], const double minv, const double maxv); + void DrawLine(gl3::GlDrawable& buff, + const double (&pts)[4][3], const double (&cv)[4], + const double minv, const double maxv); + /// Draw a 3D triangle in physical space with a central triangle removed. The /// cut is controlled by value of cut_lambda. See keys Ctrl+F3/F4. Similar to /// CutReferenceTriangle in lib/vssolution3d.cpp. diff --git a/lib/vssolution.cpp b/lib/vssolution.cpp index 7e60d869..6cc83602 100644 --- a/lib/vssolution.cpp +++ b/lib/vssolution.cpp @@ -532,7 +532,16 @@ void VisualizationSceneSolution::ToggleDrawElems() drawelems = (drawelems + 6) % 7; - cout << "Surface elements mode : " << modes[drawelems] << endl; + const int dim = mesh->Dimension(); + if (dim == 1) + { + cout << "Element mode : " << modes[drawelems] << endl; + } + else + { + cout << "Surface elements mode : " << modes[drawelems] << endl; + } + if (drawelems < 2) { extra_caption.clear(); @@ -607,8 +616,9 @@ void VisualizationSceneSolution::NewMeshAndSolution( void VisualizationSceneSolution::GetRefinedDetJ( int i, const IntegrationRule &ir, Vector &vals, DenseMatrix &tr) { - int geom = mesh->GetElementBaseGeometry(i); + const int geom = mesh->GetElementBaseGeometry(i); ElementTransformation *T = mesh->GetElementTransformation(i); + const int dim = T->GetDimension(); double Jd[4]; DenseMatrix J(Jd, 2, 2); @@ -625,15 +635,38 @@ void VisualizationSceneSolution::GetRefinedDetJ( } else if (drawelems >= 4) { - vals(j) = J.Det(); - // if (vals(j) >= 0.0) - // vals(j) = sqrt(vals(j)); - // else - // vals(j) = -sqrt(-vals(j)); + if (dim == 2) + { + vals(j) = J.Det(); + } + else if (dim == 1) + { + // Compute det(J J^T) + auto Jdl = J.Data(); + vals(j) = sqrt(Jdl[1]*Jdl[1] + Jdl[2]*Jdl[2]); // Orientation information lost. + } + else + { + std::cout << "unsupported dim " << dim << " in GetRefinedDetJ." << std::endl; + } } else { - vals(j) = J.CalcSingularvalue(0)/J.CalcSingularvalue(1); + if (dim == 2) + { + vals(j) = J.CalcSingularvalue(0)/J.CalcSingularvalue(1); + } + else if (dim == 1) + { + // Compute det(J J^T) + auto Jdl = J.Data(); + vals(j) = sqrt(Jdl[1]*Jdl[1] + Jdl[2]*Jdl[2]); // Orientation information lost. + } + else + { + std::cout << "unsupported dim " << dim << " in GetRefinedDetJ." << std::endl; + } + if (drawelems == 2) { vals(j) = vals(j) + 1.0/vals(j); @@ -687,17 +720,23 @@ int VisualizationSceneSolution::GetRefinedValuesAndNormals( { int have_normals = 0; + const int dim = mesh->Dimension(); + if (drawelems < 2) { - rsol->GetGradients(i, ir, tr); - normals.SetSize(3, tr.Width()); - for (int j = 0; j < tr.Width(); j++) + // In 1D we do not have well-defined normals. + if (dim > 1) { - normals(0, j) = -tr(0, j); - normals(1, j) = -tr(1, j); - normals(2, j) = 1.; + rsol->GetGradients(i, ir, tr); + normals.SetSize(3, tr.Width()); + for (int j = 0; j < tr.Width(); j++) + { + normals(0, j) = -tr(0, j); + normals(1, j) = -tr(1, j); + normals(2, j) = 1.; + } + have_normals = 1; } - have_normals = 1; rsol->GetValues(i, ir, vals, tr); } else @@ -1200,34 +1239,37 @@ void VisualizationSceneSolution::PrepareWithNormals() void VisualizationSceneSolution::PrepareFlat() { - int i, j; disp_buf.clear(); - int ne = mesh -> GetNE(); + const int ne = mesh -> GetNE(); DenseMatrix pointmat; Array vertices; double pts[4][3], col[4]; - for (i = 0; i < ne; i++) + for (int i = 0; i < ne; i++) { if (!el_attr_to_show[mesh->GetAttribute(i)-1]) { continue; } mesh->GetPointMatrix (i, pointmat); mesh->GetElementVertices (i, vertices); - for (j = 0; j < pointmat.Width(); j++) + for (int j = 0; j < pointmat.Width(); j++) { pts[j][0] = pointmat(0, j); pts[j][1] = pointmat(1, j); pts[j][2] = col[j] = LogVal((*sol)(vertices[j])); } - if (j == 3) + if (pointmat.Width() == 3) { DrawTriangle(disp_buf, pts, col, minv, maxv); } - else + else if (pointmat.Width() == 4) { DrawQuad(disp_buf, pts, col, minv, maxv); } + else if (pointmat.Width() == 2) + { + DrawLine(disp_buf, pts, col, minv, maxv); + } } updated_bufs.emplace_back(&disp_buf); } @@ -1241,6 +1283,8 @@ const int split_quads = 1; void VisualizationSceneSolution::PrepareFlat2() { + Array vertices; + int i, j, k; disp_buf.clear(); int ne = mesh -> GetNE(); @@ -1343,6 +1387,8 @@ void VisualizationSceneSolution::PrepareFlat2() void VisualizationSceneSolution::Prepare() { + const int dim = mesh->Dimension(); + palette.SetUseLogscale(0); switch (shading) @@ -1354,7 +1400,7 @@ void VisualizationSceneSolution::Prepare() PrepareFlat2(); return; default: - if (v_normals) + if (v_normals && dim > 1) { PrepareWithNormals(); return; @@ -1362,8 +1408,6 @@ void VisualizationSceneSolution::Prepare() break; } - int i, j; - disp_buf.clear(); gl3::GlBuilder poly = disp_buf.createBuilder(); int ne = mesh -> GetNE(); @@ -1376,47 +1420,56 @@ void VisualizationSceneSolution::Prepare() Vector ny(nv); Vector nz(nv); + // For triangles and quads: Compute the normal. + // For segments: Fill render buffer. for (int d = 0; d < mesh -> attributes.Size(); d++) { - if (!el_attr_to_show[mesh -> attributes[d]-1]) { continue; } nx = 0.; ny = 0.; nz = 0.; - for (i = 0; i < ne; i++) + // Compute normals + for (int i = 0; i < ne && (dim > 1); i++) + { if (mesh -> GetAttribute(i) == mesh -> attributes[d]) { mesh->GetPointMatrix (i, pointmat); mesh->GetElementVertices (i, vertices); - for (j = 0; j < pointmat.Size(); j++) + for (int j = 0; j < pointmat.Size(); j++) { p[j][0] = pointmat(0, j); p[j][1] = pointmat(1, j); p[j][2] = LogVal((*sol)(vertices[j])); } + int normal_state; if (pointmat.Width() == 3) { - j = Compute3DUnitNormal(p[0], p[1], p[2], nor); + normal_state = Compute3DUnitNormal(p[0], p[1], p[2], nor); } else { - j = Compute3DUnitNormal(p[0], p[1], p[2], p[3], nor); + normal_state = Compute3DUnitNormal(p[0], p[1], p[2], p[3], nor); } - if (j == 0) - for (j = 0; j < pointmat.Size(); j++) + if (normal_state == 0) // Non-degenerate normal + { + for (int j = 0; j < pointmat.Size(); j++) { nx(vertices[j]) += nor[0]; ny(vertices[j]) += nor[1]; nz(vertices[j]) += nor[2]; } + } } + } - for (i = 0; i < ne; i++) + // Fill buffers for triangles and quads. We skip this portion for + // segments, because the buffers are already filled. + for (int i = 0; i < ne; i++) { if (mesh -> GetAttribute(i) == mesh -> attributes[d]) { @@ -1429,6 +1482,9 @@ void VisualizationSceneSolution::Prepare() case Element::QUADRILATERAL: shape = GL_QUADS; break; + case Element::SEGMENT: + shape = GL_LINES; + break; default: MFEM_ABORT("Invalid 2D element type"); break; @@ -1437,7 +1493,7 @@ void VisualizationSceneSolution::Prepare() mesh->GetPointMatrix (i, pointmat); mesh->GetElementVertices (i, vertices); - for (j = 0; j < pointmat.Size(); j++) + for (int j = 0; j < pointmat.Size(); j++) { double z = LogVal((*sol)(vertices[j])); MySetColor(poly, z, minv, maxv); @@ -1486,6 +1542,12 @@ void VisualizationSceneSolution::DrawLevelCurves( gl3::GlBuilder& builder, Array &RG, DenseMatrix &pointmat, Vector &values, int sides, Array &lvl, int flat) { + const int dim = mesh->Dimension(); + if (dim == 1) // Unimplemented. + { + return; + } + double point[4][4]; // double zc = 0.5*(z[0]+z[1]); double zc = bb.z[1]; @@ -1585,7 +1647,8 @@ void VisualizationSceneSolution::PrepareLevelCurves2() void VisualizationSceneSolution::PrepareLines() { - if (shading == 2) + if (shading == 2 && + mesh->Dimension() > 1) // PrepareLines3 does not make sense for 1d meshes. { // PrepareLines2(); PrepareLines3(); @@ -1606,10 +1669,16 @@ void VisualizationSceneSolution::PrepareLines() lb.glBegin(GL_LINE_LOOP); mesh->GetPointMatrix (i, pointmat); mesh->GetElementVertices (i, vertices); - for (j = 0; j < pointmat.Size(); j++) - lb.glVertex3d(pointmat(0, j), pointmat(1, j), - LogVal((*sol)(vertices[j]))); + { + // 1D meshes get rendered flat + double z = GetMinV(); + if (mesh->Dimension() > 1) // In 1D we just put the mesh below the solution + { + z = LogVal((*sol)(vertices[j])); + } + lb.glVertex3d(pointmat(0, j), pointmat(1, j), z); + } lb.glEnd(); } @@ -1817,6 +1886,9 @@ void VisualizationSceneSolution::PrepareVertexNumbering2() void VisualizationSceneSolution::PrepareEdgeNumbering() { + // 1D meshes do not have edges. + if (mesh->Dimension() == 1) { return; } + f_nums_buf.clear(); DenseMatrix p; @@ -2012,7 +2084,6 @@ void VisualizationSceneSolution::PrepareLines3() line_buf.clear(); gl3::GlBuilder lb = line_buf.createBuilder(); - for (i = 0; i < ne; i++) { if (!el_attr_to_show[mesh->GetAttribute(i)-1]) { continue; } @@ -2063,6 +2134,12 @@ void VisualizationSceneSolution::UpdateValueRange(bool prepare) void VisualizationSceneSolution::PrepareBoundary() { + const int dim = mesh->Dimension(); + if (dim == 1) // Unimplemented. + { + return; + } + int i, j, ne = mesh->GetNBE(); Array vertices; DenseMatrix pointmat; diff --git a/lib/vssolution3d.cpp b/lib/vssolution3d.cpp index 7202a323..47c543ec 100644 --- a/lib/vssolution3d.cpp +++ b/lib/vssolution3d.cpp @@ -1486,7 +1486,7 @@ void VisualizationSceneSolution3d::PrepareFlat() DrawTriangle(disp_buf, p, c, minv, maxv); } } - else + else if (j == 4) { if (cut_lambda > 0) { @@ -1497,7 +1497,16 @@ void VisualizationSceneSolution3d::PrepareFlat() DrawQuad(disp_buf, p, c, minv, maxv); } } + else if (j == 2) + { + DrawLine(disp_buf, p, c, minv, maxv); + } + else + { + mfem_error("VisualizationSceneSolution3d::PrepareFlat() :Unknown geometry."); + } } + updated_bufs.emplace_back(&disp_buf); } @@ -1618,11 +1627,12 @@ void CutReferenceElements(int TimesToRefine, double lambda) void VisualizationSceneSolution3d::PrepareFlat2() { - int fn, fo, di, have_normals; + const int dim = mesh->Dimension(); + + int fn, fo, di, have_normals = 0; double bbox_diam, vmin, vmax; disp_buf.clear(); - int dim = mesh->Dimension(); int nbe = (dim == 3) ? mesh->GetNBE() : mesh->GetNE(); DenseMatrix pointmat, normals; Vector values, normal; @@ -1646,6 +1656,10 @@ void VisualizationSceneSolution3d::PrepareFlat2() case Element::TRIANGLE: sides = 3; break; + // TODO: can we improve this? It is quite a hack because sides + case Element::SEGMENT: + sides = 0; + break; case Element::QUADRILATERAL: default: @@ -1716,30 +1730,35 @@ void VisualizationSceneSolution3d::PrepareFlat2() CutReferenceElements(TimesToRefine, cut_lambda); cut_updated = true; } - const IntegrationRule &ir = (cut_lambda > 0) ? + const IntegrationRule &ir = (cut_lambda > 0 && dim > 1) ? ((sides == 3) ? cut_TriPts : cut_QuadPts) : RefG->RefPts; GridF->GetValues(i, ir, values, pointmat); - normals.SetSize(3, values.Size()); mesh->GetElementTransformation(i, &T); - for (int j = 0; j < values.Size(); j++) - { - T.SetIntPoint(&ir.IntPoint(j)); - const DenseMatrix &J = T.Jacobian(); - normals.GetColumnReference(j, normal); - CalcOrtho(J, normal); - normal /= normal.Norml2(); - } - have_normals = 1; di = 0; ShrinkPoints(pointmat, i, 0, 0); + + // Compute normals. Skip in 1D. + if (dim > 1) + { + normals.SetSize(3, values.Size()); + for (int j = 0; j < values.Size(); j++) + { + T.SetIntPoint(&ir.IntPoint(j)); + const DenseMatrix &J = T.Jacobian(); + normals.GetColumnReference(j, normal); + CalcOrtho(J, normal); + normal /= normal.Norml2(); + } + have_normals = 1; + } } vmin = fmin(vmin, values.Min()); vmax = fmax(vmax, values.Max()); // compute an average normal direction for the current face - if (sc != 0.0) + if (sc != 0.0 && have_normals) { for (int j = 0; j < 3; j++) { @@ -1763,7 +1782,7 @@ void VisualizationSceneSolution3d::PrepareFlat2() have_normals = 0; } - have_normals = have_normals ? 2 : 0; + have_normals = have_normals && (dim > 1) ? 2 : 0; if (di) { have_normals = -1 - have_normals; @@ -1776,6 +1795,7 @@ void VisualizationSceneSolution3d::PrepareFlat2() ((sides == 3) ? cut_TriGeoms : cut_QuadGeoms) : RefG->RefGeoms; int psides = (cut_lambda > 0) ? 4 : sides; + if (dim == 1) { psides = 2; } // Hack to trigger line rendering. DrawPatch(disp_buf, pointmat, values, normals, psides, RefGeoms, minv, maxv, have_normals); } @@ -1869,12 +1889,14 @@ void VisualizationSceneSolution3d::Prepare() { mesh->GetElementVertices(elem[i], vertices); } - for (j = 0; j < vertices.Size(); j++) + for (j = 0; j < vertices.Size() && (dim > 1); j++) { nx(vertices[j]) = ny(vertices[j]) = nz(vertices[j]) = 0.; } } - for (int i = 0; i < nelem; i++) + + // Compute normals. Skip in 1D. + for (int i = 0; i < nelem && (dim > 1); i++) { if (dim == 3) { @@ -1888,18 +1910,24 @@ void VisualizationSceneSolution3d::Prepare() } if (pointmat.Width() == 3) + { j = Compute3DUnitNormal(&pointmat(0,0), &pointmat(0,1), &pointmat(0,2), nor); + } else + { j = Compute3DUnitNormal(&pointmat(0,0), &pointmat(0,1), &pointmat(0,2), &pointmat(0,3), nor); + } if (j == 0) + { for (j = 0; j < pointmat.Size(); j++) { nx(vertices[j]) += nor[0]; ny(vertices[j]) += nor[1]; nz(vertices[j]) += nor[2]; } + } } for (int i = 0; i < nelem; i++) @@ -1936,10 +1964,14 @@ void VisualizationSceneSolution3d::Prepare() case Element::QUADRILATERAL: elemType = GL_QUADS; break; + case Element::SEGMENT: + elemType = GL_LINES; + break; default: MFEM_ABORT("Invalid boundary element type"); break; } + poly.glBegin(elemType); if (dim == 3) { @@ -1953,12 +1985,13 @@ void VisualizationSceneSolution3d::Prepare() for (j = 0; j < pointmat.Size(); j++) { MySetColor(poly, (*sol)(vertices[j]), minv, maxv); - poly.glNormal3d(nx(vertices[j]), ny(vertices[j]), nz(vertices[j])); + if (dim > 1) { poly.glNormal3d(nx(vertices[j]), ny(vertices[j]), nz(vertices[j])); } poly.glVertex3dv(&pointmat(0, j)); } poly.glEnd(); } } + updated_bufs.emplace_back(&disp_buf); } diff --git a/lib/vsvector3d.cpp b/lib/vsvector3d.cpp index 39981e38..e8986415 100644 --- a/lib/vsvector3d.cpp +++ b/lib/vsvector3d.cpp @@ -569,7 +569,7 @@ void VisualizationSceneVector3d::PrepareFlat() DrawTriangle(disp_buf, p, c, minv, maxv); } } - else + else if (j == 4) { if (cut_lambda > 0) { @@ -580,13 +580,21 @@ void VisualizationSceneVector3d::PrepareFlat() DrawQuad(disp_buf, p, c, minv, maxv); } } + else if (j == 2) + { + DrawLine(disp_buf, p, c, minv, maxv); + } + else + { + mfem_error("VisualizationSceneVector3d::PrepareFlat() :Unknown geometry."); + } } updated_bufs.emplace_back(&disp_buf); } void VisualizationSceneVector3d::PrepareFlat2() { - int fn, fo, di = 0, have_normals; + int fn, fo, di = 0, have_normals = 0; double bbox_diam, vmin, vmax; int dim = mesh->Dimension(); int ne = (dim == 3) ? mesh->GetNBE() : mesh->GetNE(); @@ -616,6 +624,10 @@ void VisualizationSceneVector3d::PrepareFlat2() case Element::TRIANGLE: sides = 3; break; + // TODO: can we improve this? It is quite a hack because sides + case Element::SEGMENT: + sides = 0; + break; case Element::QUADRILATERAL: default: @@ -685,7 +697,7 @@ void VisualizationSceneVector3d::PrepareFlat2() } ShrinkPoints(pointmat, i, fn, di); } - else // dim == 2 + else // dim < 3 { RefG = GLVisGeometryRefiner.Refine(mesh->GetElementBaseGeometry(i), TimesToRefine); @@ -695,7 +707,7 @@ void VisualizationSceneVector3d::PrepareFlat2() CutReferenceElements(TimesToRefine, cut_lambda); cut_updated = true; } - IntegrationRule &RefPts = (cut_lambda > 0) ? + IntegrationRule &RefPts = (cut_lambda > 0 && dim > 1) ? ((sides == 3) ? cut_TriPts : cut_QuadPts) : RefG->RefPts; GridF->GetValues(i, RefPts, values, pointmat); @@ -707,20 +719,24 @@ void VisualizationSceneVector3d::PrepareFlat2() } else { - const IntegrationRule &ir = (cut_lambda > 0) ? - ((sides == 3) ? cut_TriPts : cut_QuadPts) : - RefG->RefPts; - normals.SetSize(3, values.Size()); - mesh->GetElementTransformation(i, &T); - for (int j = 0; j < values.Size(); j++) + // Compute normals. Skip in 1D. + if (dim > 1) { - T.SetIntPoint(&ir.IntPoint(j)); - normals.GetColumnReference(j, normal); - CalcOrtho(T.Jacobian(), normal); - normal /= normal.Norml2(); + const IntegrationRule &ir = (cut_lambda > 0 && dim > 1) ? + ((sides == 3) ? cut_TriPts : cut_QuadPts) : + RefG->RefPts; + normals.SetSize(3, values.Size()); + mesh->GetElementTransformation(i, &T); + for (int j = 0; j < values.Size(); j++) + { + T.SetIntPoint(&ir.IntPoint(j)); + normals.GetColumnReference(j, normal); + CalcOrtho(T.Jacobian(), normal); + normal /= normal.Norml2(); + } + have_normals = 1; + di = 0; } - have_normals = 1; - di = 0; } ShrinkPoints(pointmat, i, 0, 0); } @@ -752,16 +768,17 @@ void VisualizationSceneVector3d::PrepareFlat2() have_normals = 0; } - have_normals = have_normals ? 2 : 0; + have_normals = have_normals && (dim > 1) ? 2 : 0; if (di) { have_normals = -1 - have_normals; } - Array &RefGeoms = (cut_lambda > 0) ? + Array &RefGeoms = (cut_lambda > 0 && dim > 1) ? ((sides == 3) ? cut_TriGeoms : cut_QuadGeoms) : RefG->RefGeoms; int psides = (cut_lambda > 0) ? 4 : sides; + if (dim == 1) { psides = 2; } // Hack to trigger line rendering. DrawPatch(disp_buf, pointmat, values, normals, psides, RefGeoms, minv, maxv, have_normals); } @@ -835,18 +852,24 @@ void VisualizationSceneVector3d::Prepare() } if (pointmat.Width() == 3) + { j = Compute3DUnitNormal(&pointmat(0,0), &pointmat(0,1), &pointmat(0,2), nor); + } else + { j = Compute3DUnitNormal(&pointmat(0,0), &pointmat(0,1), &pointmat(0,2), &pointmat(0,3), nor); + } if (j == 0) + { for (j = 0; j < pointmat.Size(); j++) { nx(vertices[j]) += nor[0]; ny(vertices[j]) += nor[1]; nz(vertices[j]) += nor[2]; } + } } for (i = 0; i < ne; i++) @@ -862,10 +885,15 @@ void VisualizationSceneVector3d::Prepare() case Element::TRIANGLE: draw.glBegin (GL_TRIANGLES); break; - case Element::QUADRILATERAL: draw.glBegin (GL_QUADS); break; + case Element::SEGMENT: + draw.glBegin(GL_LINES); + break; + default: + MFEM_ABORT("Invalid boundary element type"); + break; } if (dim == 3) { diff --git a/makefile b/makefile index 131efc11..7a18c3a0 100644 --- a/makefile +++ b/makefile @@ -370,7 +370,7 @@ print-%: $(info ) @true -ASTYLE_BIN = astyle +ASTYLE_BIN ?= astyle ASTYLE = $(ASTYLE_BIN) --options=$(MFEM_DIR)/config/mfem.astylerc ALL_FILES = ./glvis.cpp $(ALL_SOURCE_FILES) $(HEADER_FILES) EXT_FILES = lib/gl2ps.c lib/gl2ps.h