From ff1a6346d7586ef49ef212a7e7149b3d424179e6 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Tue, 29 Dec 2020 23:05:44 +0100 Subject: [PATCH 1/8] Turned off surface normals for integral scalar elements. --- lib/vssolution.cpp | 2 +- lib/vssolution3d.cpp | 14 +++++++++----- 2 files changed, 10 insertions(+), 6 deletions(-) diff --git a/lib/vssolution.cpp b/lib/vssolution.cpp index 0fc416a4..0f154e6c 100644 --- a/lib/vssolution.cpp +++ b/lib/vssolution.cpp @@ -725,7 +725,7 @@ int VisualizationSceneSolution::GetRefinedValuesAndNormals( if (drawelems < 2) { // In 1D we do not have well-defined normals. - if (dim > 1) + if (dim > 1 && rsol->FESpace()->GetFE(i)->GetMapType() == FiniteElement::VALUE) { rsol->GetGradients(i, ir, tr); normals.SetSize(3, tr.Width()); diff --git a/lib/vssolution3d.cpp b/lib/vssolution3d.cpp index 47c543ec..dd8bd05b 100644 --- a/lib/vssolution3d.cpp +++ b/lib/vssolution3d.cpp @@ -3915,11 +3915,7 @@ void VisualizationSceneSolution3d::PrepareLevelSurf() { RefinedGeometry *RefG; #define GLVIS_SMOOTH_LEVELSURF_NORMALS -#ifdef GLVIS_SMOOTH_LEVELSURF_NORMALS - const DenseMatrix *gp = &grad; -#else const DenseMatrix *gp = NULL; -#endif for (int ie = 0; ie < mesh->GetNE(); ie++) { @@ -3928,7 +3924,15 @@ void VisualizationSceneSolution3d::PrepareLevelSurf() RefG = GLVisGeometryRefiner.Refine(geom, TimesToRefine); GridF->GetValues(ie, RefG->RefPts, vals, pointmat); #ifdef GLVIS_SMOOTH_LEVELSURF_NORMALS - GridF->GetGradients(ie, RefG->RefPts, grad); + if (GridF->FESpace()->GetFE(ie)->GetMapType() == FiniteElement::VALUE) + { + GridF->GetGradients(ie, RefG->RefPts, grad); + gp = &grad; + } + else + { + gp = NULL; + } #endif Array &RG = RefG->RefGeoms; From 26df46d47f75e1112ab31955be2b553003fb6a83 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Wed, 5 Jun 2024 13:45:35 -0700 Subject: [PATCH 2/8] Added gradients calculation for normals. --- lib/vssolution.cpp | 29 +++++++++++++++++++++++++++-- lib/vssolution3d.cpp | 22 +++++++++++++++++++++- 2 files changed, 48 insertions(+), 3 deletions(-) diff --git a/lib/vssolution.cpp b/lib/vssolution.cpp index 0f154e6c..f0986260 100644 --- a/lib/vssolution.cpp +++ b/lib/vssolution.cpp @@ -725,9 +725,34 @@ int VisualizationSceneSolution::GetRefinedValuesAndNormals( if (drawelems < 2) { // In 1D we do not have well-defined normals. - if (dim > 1 && rsol->FESpace()->GetFE(i)->GetMapType() == FiniteElement::VALUE) + if (dim > 1) { - rsol->GetGradients(i, ir, tr); + if (rsol->FESpace()->GetFE(i)->GetMapType() == FiniteElement::VALUE) + { + rsol->GetGradients(i, ir, tr); + } + else + { + FiniteElementSpace *fes = rsol->FESpace(); + const FiniteElement *fe = fes->GetFE(i); + ElementTransformation *Trans = fes->GetElementTransformation(i); + DenseMatrix dshape(fe->GetDof(), fe->GetDim()); + Vector lval, gh(fe->GetDim()), gcol; + + rsol->GetElementDofValues(i, lval); + tr.SetSize(fe->GetDim(), ir.GetNPoints()); + for (int q = 0; q < ir.GetNPoints(); q++) + { + const IntegrationPoint &ip = ir.IntPoint(q); + fe->CalcDShape(ip, dshape); + dshape.MultTranspose(lval, gh); + Trans->SetIntPoint(&ip); + gh /= Trans->Weight(); + tr.GetColumnReference(q, gcol); + const DenseMatrix &Jinv = Trans->InverseJacobian(); + Jinv.MultTranspose(gh, gcol); + } + } normals.SetSize(3, tr.Width()); for (int j = 0; j < tr.Width(); j++) { diff --git a/lib/vssolution3d.cpp b/lib/vssolution3d.cpp index dd8bd05b..559c4877 100644 --- a/lib/vssolution3d.cpp +++ b/lib/vssolution3d.cpp @@ -3931,7 +3931,27 @@ void VisualizationSceneSolution3d::PrepareLevelSurf() } else { - gp = NULL; + FiniteElementSpace *fes = GridF->FESpace(); + const FiniteElement *fe = fes->GetFE(ie); + ElementTransformation *Trans = fes->GetElementTransformation(ie); + const IntegrationRule &ir = RefG->RefPts; + DenseMatrix dshape(fe->GetDof(), fe->GetDim()); + Vector lval, gh(fe->GetDim()), gcol; + + GridF->GetElementDofValues(ie, lval); + grad.SetSize(fe->GetDim(), ir.GetNPoints()); + for (int q = 0; q < ir.GetNPoints(); q++) + { + const IntegrationPoint &ip = ir.IntPoint(q); + fe->CalcDShape(ip, dshape); + dshape.MultTranspose(lval, gh); + Trans->SetIntPoint(&ip); + gh /= Trans->Weight(); + grad.GetColumnReference(q, gcol); + const DenseMatrix &Jinv = Trans->InverseJacobian(); + Jinv.MultTranspose(gh, gcol); + } + gp = &grad; } #endif From fb3f9adf1c91af3698e7893727d15e1ebb33137e Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Mon, 8 Jul 2024 10:19:45 -0700 Subject: [PATCH 3/8] Changed normals calculation to approximation by value-based elements. --- lib/vssolution.cpp | 18 +++++++++++++++--- lib/vssolution3d.cpp | 18 +++++++++++++++--- 2 files changed, 30 insertions(+), 6 deletions(-) diff --git a/lib/vssolution.cpp b/lib/vssolution.cpp index f0986260..d7f61f78 100644 --- a/lib/vssolution.cpp +++ b/lib/vssolution.cpp @@ -735,11 +735,24 @@ int VisualizationSceneSolution::GetRefinedValuesAndNormals( { FiniteElementSpace *fes = rsol->FESpace(); const FiniteElement *fe = fes->GetFE(i); + const int ndof = fe->GetDof(); + const int ndim = fe->GetDim(); ElementTransformation *Trans = fes->GetElementTransformation(i); - DenseMatrix dshape(fe->GetDof(), fe->GetDim()); - Vector lval, gh(fe->GetDim()), gcol; + DenseMatrix dshape(ndof, ndim); + Vector lval, gh(ndim), gcol; rsol->GetElementDofValues(i, lval); + + // Local projection to value-based FE + const IntegrationRule &nodes = fe->GetNodes(); + for (int n = 0; n < nodes.GetNPoints(); n++) + { + const IntegrationPoint &ip = nodes.IntPoint(n); + Trans->SetIntPoint(&ip); + lval(n) /= Trans->Weight(); + } + + // Gradient calculation tr.SetSize(fe->GetDim(), ir.GetNPoints()); for (int q = 0; q < ir.GetNPoints(); q++) { @@ -747,7 +760,6 @@ int VisualizationSceneSolution::GetRefinedValuesAndNormals( fe->CalcDShape(ip, dshape); dshape.MultTranspose(lval, gh); Trans->SetIntPoint(&ip); - gh /= Trans->Weight(); tr.GetColumnReference(q, gcol); const DenseMatrix &Jinv = Trans->InverseJacobian(); Jinv.MultTranspose(gh, gcol); diff --git a/lib/vssolution3d.cpp b/lib/vssolution3d.cpp index 559c4877..ff21ebd9 100644 --- a/lib/vssolution3d.cpp +++ b/lib/vssolution3d.cpp @@ -3933,12 +3933,25 @@ void VisualizationSceneSolution3d::PrepareLevelSurf() { FiniteElementSpace *fes = GridF->FESpace(); const FiniteElement *fe = fes->GetFE(ie); + const int ndof = fe->GetDof(); + const int ndim = fe->GetDim(); ElementTransformation *Trans = fes->GetElementTransformation(ie); const IntegrationRule &ir = RefG->RefPts; - DenseMatrix dshape(fe->GetDof(), fe->GetDim()); - Vector lval, gh(fe->GetDim()), gcol; + DenseMatrix dshape(ndof, ndim); + Vector lval, gh(ndim), gcol; GridF->GetElementDofValues(ie, lval); + + // Local projection to value-based FE + const IntegrationRule &nodes = fe->GetNodes(); + for (int n = 0; n < nodes.GetNPoints(); n++) + { + const IntegrationPoint &ip = nodes.IntPoint(n); + Trans->SetIntPoint(&ip); + lval(n) /= Trans->Weight(); + } + + // Gradient calculation grad.SetSize(fe->GetDim(), ir.GetNPoints()); for (int q = 0; q < ir.GetNPoints(); q++) { @@ -3946,7 +3959,6 @@ void VisualizationSceneSolution3d::PrepareLevelSurf() fe->CalcDShape(ip, dshape); dshape.MultTranspose(lval, gh); Trans->SetIntPoint(&ip); - gh /= Trans->Weight(); grad.GetColumnReference(q, gcol); const DenseMatrix &Jinv = Trans->InverseJacobian(); Jinv.MultTranspose(gh, gcol); From ca3624c6c8309782b8e836d137da95f492252a84 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Mon, 8 Jul 2024 15:38:29 -0700 Subject: [PATCH 4/8] Added a note to Changelog. --- CHANGELOG | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/CHANGELOG b/CHANGELOG index e42bf778..a70fbff5 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -16,6 +16,10 @@ Version 4.2.1 (development) - Add support to visualize solutions on 1D elements embedded in 2D and 3D. +- Added support for scalar integral finite elements, where calculation of the + surface normals was not implemented and was crashing GLVis. The normals are + approximately calculated from the point-wise projected value-based FEs. + Version 4.2 released on May 23, 2022 ==================================== From f6bb46c527c5580c20920c417cb534c72310f1ae Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Tue, 9 Jul 2024 09:22:43 -0700 Subject: [PATCH 5/8] Added small comments about local projection. --- lib/vssolution.cpp | 2 +- lib/vssolution3d.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/vssolution.cpp b/lib/vssolution.cpp index d7f61f78..9f0e90fc 100644 --- a/lib/vssolution.cpp +++ b/lib/vssolution.cpp @@ -749,7 +749,7 @@ int VisualizationSceneSolution::GetRefinedValuesAndNormals( { const IntegrationPoint &ip = nodes.IntPoint(n); Trans->SetIntPoint(&ip); - lval(n) /= Trans->Weight(); + lval(n) /= Trans->Weight();//value = dof / |J| } // Gradient calculation diff --git a/lib/vssolution3d.cpp b/lib/vssolution3d.cpp index ff21ebd9..a3db1962 100644 --- a/lib/vssolution3d.cpp +++ b/lib/vssolution3d.cpp @@ -3948,7 +3948,7 @@ void VisualizationSceneSolution3d::PrepareLevelSurf() { const IntegrationPoint &ip = nodes.IntPoint(n); Trans->SetIntPoint(&ip); - lval(n) /= Trans->Weight(); + lval(n) /= Trans->Weight();//value = dof / |J| } // Gradient calculation From f540a98346b4dcabaf3e3b6a67a98757b73257ed Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Tue, 9 Jul 2024 09:29:20 -0700 Subject: [PATCH 6/8] More explicit mapping type check. --- lib/vssolution.cpp | 10 ++++++++-- lib/vssolution3d.cpp | 9 +++++++-- 2 files changed, 15 insertions(+), 4 deletions(-) diff --git a/lib/vssolution.cpp b/lib/vssolution.cpp index 9f0e90fc..944db1d6 100644 --- a/lib/vssolution.cpp +++ b/lib/vssolution.cpp @@ -727,11 +727,12 @@ int VisualizationSceneSolution::GetRefinedValuesAndNormals( // In 1D we do not have well-defined normals. if (dim > 1) { - if (rsol->FESpace()->GetFE(i)->GetMapType() == FiniteElement::VALUE) + const int map_type = rsol->FESpace()->GetFE(i)->GetMapType(); + if (map_type == FiniteElement::MapType::VALUE) { rsol->GetGradients(i, ir, tr); } - else + else (map_type == FiniteElement::MapType::INTEGRAL) { FiniteElementSpace *fes = rsol->FESpace(); const FiniteElement *fe = fes->GetFE(i); @@ -765,6 +766,11 @@ int VisualizationSceneSolution::GetRefinedValuesAndNormals( Jinv.MultTranspose(gh, gcol); } } + else + { + MFEM_ABORT("Unknown mapping type"); + } + normals.SetSize(3, tr.Width()); for (int j = 0; j < tr.Width(); j++) { diff --git a/lib/vssolution3d.cpp b/lib/vssolution3d.cpp index a3db1962..2a0daf52 100644 --- a/lib/vssolution3d.cpp +++ b/lib/vssolution3d.cpp @@ -3924,12 +3924,13 @@ void VisualizationSceneSolution3d::PrepareLevelSurf() RefG = GLVisGeometryRefiner.Refine(geom, TimesToRefine); GridF->GetValues(ie, RefG->RefPts, vals, pointmat); #ifdef GLVIS_SMOOTH_LEVELSURF_NORMALS - if (GridF->FESpace()->GetFE(ie)->GetMapType() == FiniteElement::VALUE) + const int map_type = GridF->FESpace()->GetFE(ie)->GetMapType(); + if (map_type == FiniteElement::MapType::VALUE) { GridF->GetGradients(ie, RefG->RefPts, grad); gp = &grad; } - else + else if (map_type == FiniteElement::MapType::INTEGRAL) { FiniteElementSpace *fes = GridF->FESpace(); const FiniteElement *fe = fes->GetFE(ie); @@ -3965,6 +3966,10 @@ void VisualizationSceneSolution3d::PrepareLevelSurf() } gp = &grad; } + else + { + MFEM_ABORT("Unknown mapping type"); + } #endif Array &RG = RefG->RefGeoms; From 99d236e9f5ca12cd2ae31ff7fadf7d8a1b331fa7 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Tue, 9 Jul 2024 09:34:01 -0700 Subject: [PATCH 7/8] Fixed typo in if statement. --- lib/vssolution.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/vssolution.cpp b/lib/vssolution.cpp index 944db1d6..91e853eb 100644 --- a/lib/vssolution.cpp +++ b/lib/vssolution.cpp @@ -732,7 +732,7 @@ int VisualizationSceneSolution::GetRefinedValuesAndNormals( { rsol->GetGradients(i, ir, tr); } - else (map_type == FiniteElement::MapType::INTEGRAL) + else if (map_type == FiniteElement::MapType::INTEGRAL) { FiniteElementSpace *fes = rsol->FESpace(); const FiniteElement *fe = fes->GetFE(i); From 380c2e9f30e2e04adce94d07b4af970c0b9b9fcc Mon Sep 17 00:00:00 2001 From: Tzanio Kolev Date: Wed, 17 Jul 2024 10:40:33 -0700 Subject: [PATCH 8/8] minor --- CHANGELOG | 2 +- tests/data | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/CHANGELOG b/CHANGELOG index 16adf868..7d6fab7a 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -18,7 +18,7 @@ Version 4.2.1 (development) - Added support for scalar integral finite elements, where calculation of the surface normals was not implemented and was crashing GLVis. The normals are - approximately calculated from the point-wise projected value-based FEs. + approximately calculated from the point-wise projected value-based elements. - Added two new modes for visualization of vector fields in 2D, placing the arrows above the plotted surface and using a single color. diff --git a/tests/data b/tests/data index 81de802a..53fbee22 160000 --- a/tests/data +++ b/tests/data @@ -1 +1 @@ -Subproject commit 81de802ad63ba41cd1fcd3825fbb3d2471538ace +Subproject commit 53fbee222cceb0581122cb72518bfc1b85766b4e