From 84a101a745381ec544d8e65ae0605a60a1b8e053 Mon Sep 17 00:00:00 2001 From: Will Pazner Date: Mon, 14 Aug 2023 17:28:34 +0900 Subject: [PATCH 1/5] Add 1D meshes embedded in 2D and 3D --- data/diag-segment-2d.mesh | 37 +++++++++++++++++++++++++++++++++++++ data/diag-segment-3d.mesh | 37 +++++++++++++++++++++++++++++++++++++ 2 files changed, 74 insertions(+) create mode 100644 data/diag-segment-2d.mesh create mode 100644 data/diag-segment-3d.mesh diff --git a/data/diag-segment-2d.mesh b/data/diag-segment-2d.mesh new file mode 100644 index 00000000000..793d883d85e --- /dev/null +++ b/data/diag-segment-2d.mesh @@ -0,0 +1,37 @@ +MFEM mesh v1.0 + +# +# MFEM Geometry Types (see mesh/geom.hpp): +# +# POINT = 0 +# SEGMENT = 1 +# TRIANGLE = 2 +# SQUARE = 3 +# TETRAHEDRON = 4 +# CUBE = 5 +# PRISM = 6 +# + +dimension +1 + +elements +4 +1 1 0 1 +1 1 1 2 +1 1 2 3 +1 1 3 4 + +boundary +2 +1 0 0 +2 0 4 + +vertices +5 +2 +0 0 +0.25 0.25 +0.50 0.50 +0.75 0.75 +1 1 diff --git a/data/diag-segment-3d.mesh b/data/diag-segment-3d.mesh new file mode 100644 index 00000000000..593c609986d --- /dev/null +++ b/data/diag-segment-3d.mesh @@ -0,0 +1,37 @@ +MFEM mesh v1.0 + +# +# MFEM Geometry Types (see mesh/geom.hpp): +# +# POINT = 0 +# SEGMENT = 1 +# TRIANGLE = 2 +# SQUARE = 3 +# TETRAHEDRON = 4 +# CUBE = 5 +# PRISM = 6 +# + +dimension +1 + +elements +4 +1 1 0 1 +1 1 1 2 +1 1 2 3 +1 1 3 4 + +boundary +2 +1 0 0 +2 0 4 + +vertices +5 +3 +0 0 0 +0.25 0.25 0.25 +0.50 0.50 0.50 +0.75 0.75 0.75 +1 1 1 From a0f5bd4e444c22dd44e9276db757d20c5fb73aac Mon Sep 17 00:00:00 2001 From: Will Pazner Date: Mon, 14 Aug 2023 17:31:16 +0900 Subject: [PATCH 2/5] Support sdim != dim (surface meshes) in QuadratureInterpolator::PhysDerivatives --- fem/qinterp/grad.hpp | 129 ++++++++++++++++++++++------- fem/qinterp/grad_by_nodes.cpp | 4 +- fem/qinterp/grad_by_vdim.cpp | 4 +- fem/qinterp/grad_phys_by_nodes.cpp | 30 +++---- fem/qinterp/grad_phys_by_vdim.cpp | 20 ++--- 5 files changed, 126 insertions(+), 61 deletions(-) diff --git a/fem/qinterp/grad.hpp b/fem/qinterp/grad.hpp index 7c9f6db503d..61228b8aace 100644 --- a/fem/qinterp/grad.hpp +++ b/fem/qinterp/grad.hpp @@ -33,16 +33,17 @@ static void Derivatives1D(const int NE, const double *j_, const double *x_, double *y_, + const int sdim, const int vdim, const int d1d, const int q1d) { const auto g = Reshape(g_, q1d, d1d); - const auto j = Reshape(j_, q1d, NE); + const auto j = Reshape(j_, q1d, sdim, NE); const auto x = Reshape(x_, d1d, vdim, NE); auto y = Q_LAYOUT == QVectorLayout::byNODES ? - Reshape(y_, q1d, vdim, NE): - Reshape(y_, vdim, q1d, NE); + Reshape(y_, q1d, vdim, sdim, NE): + Reshape(y_, vdim, sdim, q1d, NE); mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e) { @@ -50,14 +51,40 @@ static void Derivatives1D(const int NE, { for (int q = 0; q < q1d; q++) { - double u = 0.0; + double du[3] = {0.0, 0.0, 0.0}; for (int d = 0; d < d1d; d++) { - u += g(q, d) * x(d, c, e); + du[0] += g(q, d) * x(d, c, e); + } + if (GRAD_PHYS) + { + if (sdim == 1) { du[0] /= j(q, 0, e); } + else if (sdim == 2) + { + const double Jloc[2] = {j(q,0,e), j(q,1,e)}; + const double t = 1.0 / (Jloc[0]*Jloc[0] + Jloc[1]*Jloc[1]); + const double U = t*Jloc[0]*du[0]; + const double V = t*Jloc[1]*du[0]; + du[0] = U; + du[1] = V; + } + else // sdim == 3 + { + const double Jloc[3] = {j(q,0,e), j(q,1,e), j(q,2,e)}; + const double t = 1.0/(Jloc[0]*Jloc[0] + Jloc[1]*Jloc[1] + Jloc[2]*Jloc[2]); + const double U = t*Jloc[0]*du[0]; + const double V = t*Jloc[1]*du[0]; + const double W = t*Jloc[2]*du[0]; + du[0] = U; + du[1] = V; + du[2] = W; + } + } + for (int d = 0; d < sdim; ++d) + { + if (Q_LAYOUT == QVectorLayout::byVDIM) { y(c, d, q, e) = du[d]; } + if (Q_LAYOUT == QVectorLayout::byNODES) { y(q, c, d, e) = du[d]; } } - if (GRAD_PHYS) { u /= j(q, e); } - if (Q_LAYOUT == QVectorLayout::byVDIM) { y(c, q, e) = u; } - if (Q_LAYOUT == QVectorLayout::byNODES) { y(q, c, e) = u; } } } }); @@ -73,6 +100,7 @@ static void Derivatives2D(const int NE, const double *j_, const double *x_, double *y_, + const int sdim = 2, const int vdim = 0, const int d1d = 0, const int q1d = 0) @@ -80,15 +108,16 @@ static void Derivatives2D(const int NE, const int D1D = T_D1D ? T_D1D : d1d; const int Q1D = T_Q1D ? T_Q1D : q1d; const int VDIM = T_VDIM ? T_VDIM : vdim; + const int SDIM = GRAD_PHYS ? sdim : 2; static constexpr int NBZ = T_NBZ ? T_NBZ : 1; const auto b = Reshape(b_, Q1D, D1D); const auto g = Reshape(g_, Q1D, D1D); - const auto j = Reshape(j_, Q1D, Q1D, 2, 2, NE); + const auto j = Reshape(j_, Q1D, Q1D, SDIM, 2, NE); const auto x = Reshape(x_, D1D, D1D, VDIM, NE); auto y = Q_LAYOUT == QVectorLayout:: byNODES ? - Reshape(y_, Q1D, Q1D, VDIM, 2, NE): - Reshape(y_, VDIM, 2, Q1D, Q1D, NE); + Reshape(y_, Q1D, Q1D, VDIM, SDIM, NE): + Reshape(y_, VDIM, SDIM, Q1D, Q1D, NE); mfem::forall_2D_batch(NE, Q1D, Q1D, NBZ, [=] MFEM_HOST_DEVICE (int e) { @@ -135,34 +164,70 @@ static void Derivatives2D(const int NE, { MFEM_FOREACH_THREAD(qx,x,Q1D) { - double u = 0.0; - double v = 0.0; + double du[3] = {0.0, 0.0, 0.0}; for (int dy = 0; dy < D1D; ++dy) { - u += DQ1(dy,qx) * B(dy,qy); - v += DQ0(dy,qx) * G(dy,qy); + du[0] += DQ1(dy,qx) * B(dy,qy); + du[1] += DQ0(dy,qx) * G(dy,qy); } if (GRAD_PHYS) { - double Jloc[4], Jinv[4]; - Jloc[0] = j(qx,qy,0,0,e); - Jloc[1] = j(qx,qy,1,0,e); - Jloc[2] = j(qx,qy,0,1,e); - Jloc[3] = j(qx,qy,1,1,e); - kernels::CalcInverse<2>(Jloc, Jinv); - const double U = Jinv[0]*u + Jinv[1]*v; - const double V = Jinv[2]*u + Jinv[3]*v; - u = U; v = V; - } - if (Q_LAYOUT == QVectorLayout::byVDIM) - { - y(c,0,qx,qy,e) = u; - y(c,1,qx,qy,e) = v; + if (SDIM == 2) + { + double Jloc[4], Jinv[4]; + Jloc[0] = j(qx,qy,0,0,e); + Jloc[1] = j(qx,qy,1,0,e); + Jloc[2] = j(qx,qy,0,1,e); + Jloc[3] = j(qx,qy,1,1,e); + kernels::CalcInverse<2>(Jloc, Jinv); + const double U = Jinv[0]*du[0] + Jinv[1]*du[1]; + const double V = Jinv[2]*du[0] + Jinv[3]*du[1]; + du[0] = U; + du[1] = V; + } + else + { + double Jloc[6], Jinv[6]; + Jloc[0] = j(qx,qy,0,0,e); + Jloc[1] = j(qx,qy,1,0,e); + Jloc[2] = j(qx,qy,2,0,e); + Jloc[3] = j(qx,qy,0,1,e); + Jloc[4] = j(qx,qy,1,1,e); + Jloc[5] = j(qx,qy,2,1,e); + + double ee, gg, ff; + ee = Jloc[0]*Jloc[0] + Jloc[1]*Jloc[1] + Jloc[2]*Jloc[2]; + gg = Jloc[3]*Jloc[3] + Jloc[4]*Jloc[4] + Jloc[5]*Jloc[5]; + ff = Jloc[0]*Jloc[3] + Jloc[1]*Jloc[4] + Jloc[2]*Jloc[5]; + const double t = 1.0 / (ee*gg - ff*ff); + ee *= t; gg *= t; ff *= t; + + Jinv[0] = Jloc[0]*gg - Jloc[3]*ff; + Jinv[1] = Jloc[3]*ee - Jloc[0]*ff; + Jinv[2] = Jloc[1]*gg - Jloc[4]*ff; + Jinv[3] = Jloc[4]*ee - Jloc[1]*ff; + Jinv[4] = Jloc[2]*gg - Jloc[5]*ff; + Jinv[5] = Jloc[5]*ee - Jloc[2]*ff; + + const double U = Jinv[0]*du[0] + Jinv[1]*du[1]; + const double V = Jinv[2]*du[0] + Jinv[3]*du[1]; + const double W = Jinv[4]*du[0] + Jinv[5]*du[1]; + + du[0] = U; + du[1] = V; + du[2] = W; + } } - if (Q_LAYOUT == QVectorLayout::byNODES) + for (int d = 0; d < SDIM; ++d) { - y(qx,qy,c,0,e) = u; - y(qx,qy,c,1,e) = v; + if (Q_LAYOUT == QVectorLayout::byVDIM) + { + y(c,d,qx,qy,e) = du[d]; + } + else // Q_LAYOUT == QVectorLayout::byNODES + { + y(qx,qy,c,d,e) = du[d]; + } } } } diff --git a/fem/qinterp/grad_by_nodes.cpp b/fem/qinterp/grad_by_nodes.cpp index 2b8ef92ce43..80bd6a04714 100644 --- a/fem/qinterp/grad_by_nodes.cpp +++ b/fem/qinterp/grad_by_nodes.cpp @@ -47,7 +47,7 @@ void TensorDerivatives(const int NE, if (dim == 1) { - return Derivatives1D(NE,G,J,X,Y,vdim,D1D,Q1D); + return Derivatives1D(NE,G,J,X,Y,dim,vdim,D1D,Q1D); } if (dim == 2) { @@ -83,7 +83,7 @@ void TensorDerivatives(const int NE, { MFEM_ABORT(""); } - Derivatives2D(NE,B,G,J,X,Y,vdim,D1D,Q1D); + Derivatives2D(NE,B,G,J,X,Y,dim,vdim,D1D,Q1D); return; } } diff --git a/fem/qinterp/grad_by_vdim.cpp b/fem/qinterp/grad_by_vdim.cpp index 15eb7be2e65..32902dfa9f8 100644 --- a/fem/qinterp/grad_by_vdim.cpp +++ b/fem/qinterp/grad_by_vdim.cpp @@ -47,7 +47,7 @@ void TensorDerivatives(const int NE, if (dim == 1) { - return Derivatives1D(NE,G,J,X,Y,vdim,D1D,Q1D); + return Derivatives1D(NE,G,J,X,Y,dim,vdim,D1D,Q1D); } if (dim == 2) { @@ -68,7 +68,7 @@ void TensorDerivatives(const int NE, << " are not supported!"); MFEM_VERIFY(Q1D <= MQ, "Quadrature rules with more than " << MQ << " 1D points are not supported!"); - Derivatives2D(NE,B,G,J,X,Y,vdim,D1D,Q1D); + Derivatives2D(NE,B,G,J,X,Y,dim,vdim,D1D,Q1D); return; } } diff --git a/fem/qinterp/grad_phys_by_nodes.cpp b/fem/qinterp/grad_phys_by_nodes.cpp index 9ce6a20e23f..0d1cf16a799 100644 --- a/fem/qinterp/grad_phys_by_nodes.cpp +++ b/fem/qinterp/grad_phys_by_nodes.cpp @@ -37,7 +37,7 @@ void TensorPhysDerivatives(const int NE, const int D1D = maps.ndof; const int Q1D = maps.nqpt; - MFEM_ASSERT(geom.mesh->SpaceDimension() == dim, ""); + const int sdim = geom.mesh->SpaceDimension(); const double *B = maps.B.Read(); const double *G = maps.G.Read(); @@ -52,25 +52,25 @@ void TensorPhysDerivatives(const int NE, if (dim == 1) { - return Derivatives1D(NE,G,J,X,Y,vdim,D1D,Q1D); + return Derivatives1D(NE,G,J,X,Y,sdim,vdim,D1D,Q1D); } if (dim == 2) { switch (id) { - case 0x133: return Derivatives2D(NE,B,G,J,X,Y); - case 0x134: return Derivatives2D(NE,B,G,J,X,Y); - case 0x143: return Derivatives2D(NE,B,G,J,X,Y); - case 0x144: return Derivatives2D(NE,B,G,J,X,Y); - case 0x146: return Derivatives2D(NE,B,G,J,X,Y); - case 0x158: return Derivatives2D(NE,B,G,J,X,Y); + case 0x133: return Derivatives2D(NE,B,G,J,X,Y,sdim); + case 0x134: return Derivatives2D(NE,B,G,J,X,Y,sdim); + case 0x143: return Derivatives2D(NE,B,G,J,X,Y,sdim); + case 0x144: return Derivatives2D(NE,B,G,J,X,Y,sdim); + case 0x146: return Derivatives2D(NE,B,G,J,X,Y,sdim); + case 0x158: return Derivatives2D(NE,B,G,J,X,Y,sdim); - case 0x233: return Derivatives2D(NE,B,G,J,X,Y); - case 0x234: return Derivatives2D(NE,B,G,J,X,Y); - case 0x243: return Derivatives2D(NE,B,G,J,X,Y); - case 0x244: return Derivatives2D(NE,B,G,J,X,Y); - case 0x246: return Derivatives2D(NE,B,G,J,X,Y); - case 0x258: return Derivatives2D(NE,B,G,J,X,Y); + case 0x233: return Derivatives2D(NE,B,G,J,X,Y,sdim); + case 0x234: return Derivatives2D(NE,B,G,J,X,Y,sdim); + case 0x243: return Derivatives2D(NE,B,G,J,X,Y,sdim); + case 0x244: return Derivatives2D(NE,B,G,J,X,Y,sdim); + case 0x246: return Derivatives2D(NE,B,G,J,X,Y,sdim); + case 0x258: return Derivatives2D(NE,B,G,J,X,Y,sdim); default: { constexpr int MD = MAX_D1D; @@ -79,7 +79,7 @@ void TensorPhysDerivatives(const int NE, << " are not supported!"); MFEM_VERIFY(Q1D <= MQ, "Quadrature rules with more than " << MQ << " 1D points are not supported!"); - Derivatives2D(NE,B,G,J,X,Y,vdim,D1D,Q1D); + Derivatives2D(NE,B,G,J,X,Y,sdim,vdim,D1D,Q1D); return; } } diff --git a/fem/qinterp/grad_phys_by_vdim.cpp b/fem/qinterp/grad_phys_by_vdim.cpp index 8bb9370997c..6004ae14e04 100644 --- a/fem/qinterp/grad_phys_by_vdim.cpp +++ b/fem/qinterp/grad_phys_by_vdim.cpp @@ -37,7 +37,7 @@ void TensorPhysDerivatives(const int NE, const int D1D = maps.ndof; const int Q1D = maps.nqpt; - MFEM_ASSERT(geom.mesh->SpaceDimension() == dim, ""); + const int sdim = geom.mesh->SpaceDimension(); const double *B = maps.B.Read(); const double *G = maps.G.Read(); @@ -52,20 +52,20 @@ void TensorPhysDerivatives(const int NE, if (dim == 1) { - return Derivatives1D(NE,G,J,X,Y,vdim,D1D,Q1D); + return Derivatives1D(NE,G,J,X,Y,sdim,vdim,D1D,Q1D); } if (dim == 2) { switch (id) { - case 0x134: return Derivatives2D(NE,B,G,J,X,Y); - case 0x146: return Derivatives2D(NE,B,G,J,X,Y); - case 0x158: return Derivatives2D(NE,B,G,J,X,Y); + case 0x134: return Derivatives2D(NE,B,G,J,X,Y,sdim); + case 0x146: return Derivatives2D(NE,B,G,J,X,Y,sdim); + case 0x158: return Derivatives2D(NE,B,G,J,X,Y,sdim); - case 0x233: return Derivatives2D(NE,B,G,J,X,Y); - case 0x234: return Derivatives2D(NE,B,G,J,X,Y); - case 0x246: return Derivatives2D(NE,B,G,J,X,Y); - case 0x258: return Derivatives2D(NE,B,G,J,X,Y); + case 0x233: return Derivatives2D(NE,B,G,J,X,Y,sdim); + case 0x234: return Derivatives2D(NE,B,G,J,X,Y,sdim); + case 0x246: return Derivatives2D(NE,B,G,J,X,Y,sdim); + case 0x258: return Derivatives2D(NE,B,G,J,X,Y,sdim); default: { constexpr int MD = MAX_D1D; @@ -74,7 +74,7 @@ void TensorPhysDerivatives(const int NE, << " are not supported!"); MFEM_VERIFY(Q1D <= MQ, "Quadrature rules with more than " << MQ << " 1D points are not supported!"); - Derivatives2D(NE,B,G,J,X,Y,vdim,D1D,Q1D); + Derivatives2D(NE,B,G,J,X,Y,sdim,vdim,D1D,Q1D); return; } } From 4d94ce5ac8b0b4d48c64e6f79a211cb60b0a8d56 Mon Sep 17 00:00:00 2001 From: Will Pazner Date: Mon, 14 Aug 2023 17:31:50 +0900 Subject: [PATCH 3/5] Add more QuadratureInterpolator unit tests --- tests/unit/fem/test_quadinterpolator.cpp | 112 ++++++++++++++++++++--- 1 file changed, 100 insertions(+), 12 deletions(-) diff --git a/tests/unit/fem/test_quadinterpolator.cpp b/tests/unit/fem/test_quadinterpolator.cpp index e68c2f0bfab..f9862c8a08d 100644 --- a/tests/unit/fem/test_quadinterpolator.cpp +++ b/tests/unit/fem/test_quadinterpolator.cpp @@ -236,16 +236,104 @@ static bool testQuadratureInterpolator(const int dim, return true; } -TEST_CASE("QuadratureInterpolator", - "[QuadratureInterpolator]" - "[CUDA]") +TEST_CASE("QuadratureInterpolator", "[QuadratureInterpolator][CUDA]") { - const auto dim = GENERATE(1,2,3); // dimension - const auto p = GENERATE(range(1,7)); // element order, 1 <= p < 7 - const auto q = GENERATE_COPY(p+1,p+2); // 1D quadrature points - const auto l = GENERATE(QVectorLayout::byNODES, QVectorLayout::byVDIM); - const auto nx = 3; // number of element in x - const auto ny = 3; // number of element in y - const auto nz = 3; // number of element in z - testQuadratureInterpolator(dim, p, q, l, nx, ny, nz); -} // TEST_CASE "QuadratureInterpolator" + SECTION("Tensor and non-tensor") + { + const auto dim = GENERATE(1,2,3); // dimension + const auto p = GENERATE(range(1,7)); // element order, 1 <= p < 7 + const auto q = GENERATE_COPY(p+1,p+2); // 1D quadrature points + const auto l = GENERATE(QVectorLayout::byNODES, QVectorLayout::byVDIM); + const auto nx = 3; // number of element in x + const auto ny = 3; // number of element in y + const auto nz = 3; // number of element in z + testQuadratureInterpolator(dim, p, q, l, nx, ny, nz); + } + + SECTION("Values and physical derivatives") + { + const auto mesh_fname = GENERATE( + "../../data/inline-segment.mesh", + "../../data/star.mesh", + "../../data/star-q3.mesh", + "../../data/fichera.mesh", + "../../data/fichera-q3.mesh", + "../../data/diag-segment-2d.mesh", // 1D mesh in 2D + "../../data/diag-segment-3d.mesh", // 1D mesh in 3D + "../../data/star-surf.mesh" // surface mesh + ); + const int order = GENERATE(1, 2, 3); + + Mesh mesh = Mesh::LoadFromFile(mesh_fname); + H1_FECollection fec(order); + FiniteElementSpace fes(&mesh, &fec); + QuadratureSpace qs(&mesh, 2*order); + + GridFunction gf(&fes); + gf.Randomize(1); + + const ElementDofOrdering ordering = ElementDofOrdering::LEXICOGRAPHIC; + // Use element restriction to go from L-vector to E-vector + const Operator *R = fes.GetElementRestriction(ordering); + Vector e_vec(R->Height()); + R->Mult(gf, e_vec); + + // Use quadrature interpolator to go from E-vector to Q-vector + const QuadratureInterpolator *qi = + fes.GetQuadratureInterpolator(qs); + qi->SetOutputLayout(QVectorLayout::byVDIM); + + { + QuadratureFunction qf1(qs), qf2(qs); + + const int ne = qs.GetNE(); + Vector values; + for (int iel = 0; iel < ne; ++iel) + { + qf1.GetValues(iel, values); + const IntegrationRule &ir = qs.GetIntRule(iel); + ElementTransformation &T = *qs.GetTransformation(iel); + for (int iq = 0; iq < ir.Size(); ++iq) + { + const IntegrationPoint &ip = ir[iq]; + T.SetIntPoint(&ip); + const int iq_p = qs.GetPermutedIndex(iel, iq); + values[iq_p] = gf.GetValue(T, ip); + } + } + + qi->Values(e_vec, qf2); + + qf1 -= qf2; + REQUIRE(qf1.Normlinf() == MFEM_Approx(0.0)); + } + + { + const int sdim = mesh.SpaceDimension(); + QuadratureFunction qf1(qs, sdim), qf2(qs, sdim); + + const int ne = qs.GetNE(); + DenseMatrix values; + Vector col; + for (int iel = 0; iel < ne; ++iel) + { + qf1.GetValues(iel, values); + const IntegrationRule &ir = qs.GetIntRule(iel); + ElementTransformation &T = *qs.GetTransformation(iel); + for (int iq = 0; iq < ir.Size(); ++iq) + { + const IntegrationPoint &ip = ir[iq]; + T.SetIntPoint(&ip); + const int iq_p = qs.GetPermutedIndex(iel, iq); + values.GetColumnReference(iq_p, col); + gf.GetGradient(T, col); + } + } + + qi->PhysDerivatives(e_vec, qf2); + + qf1 -= qf2; + REQUIRE(qf1.Normlinf() == MFEM_Approx(0.0)); + } + } +} From 7284dc092fbc3fb8d7f9d9a73834814c8b511d39 Mon Sep 17 00:00:00 2001 From: Will Pazner Date: Tue, 12 Sep 2023 11:54:58 -0700 Subject: [PATCH 4/5] Remove un-needed include --- linalg/densemat.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/linalg/densemat.cpp b/linalg/densemat.cpp index 44b40fa82e1..15ef7c0cfc1 100644 --- a/linalg/densemat.cpp +++ b/linalg/densemat.cpp @@ -17,7 +17,6 @@ #include "vector.hpp" #include "matrix.hpp" #include "densemat.hpp" -#include "kernels.hpp" #include "../general/forall.hpp" #include "../general/table.hpp" #include "../general/globals.hpp" From 6b5db0c503c8899930329332769f12011006ed4c Mon Sep 17 00:00:00 2001 From: Will Pazner Date: Tue, 12 Sep 2023 11:55:15 -0700 Subject: [PATCH 5/5] Add linalg kernels::CalcLeftInverse --- fem/qinterp/grad.hpp | 33 ++++++++++----------------------- linalg/densemat.cpp | 27 ++++----------------------- linalg/kernels.hpp | 39 +++++++++++++++++++++++++++++++++++++++ 3 files changed, 53 insertions(+), 46 deletions(-) diff --git a/fem/qinterp/grad.hpp b/fem/qinterp/grad.hpp index 61228b8aace..8f0a2d3d2a0 100644 --- a/fem/qinterp/grad.hpp +++ b/fem/qinterp/grad.hpp @@ -62,19 +62,21 @@ static void Derivatives1D(const int NE, else if (sdim == 2) { const double Jloc[2] = {j(q,0,e), j(q,1,e)}; - const double t = 1.0 / (Jloc[0]*Jloc[0] + Jloc[1]*Jloc[1]); - const double U = t*Jloc[0]*du[0]; - const double V = t*Jloc[1]*du[0]; + double Jinv[3]; + kernels::CalcLeftInverse<2,1>(Jloc, Jinv); + const double U = Jinv[0]*du[0]; + const double V = Jinv[1]*du[0]; du[0] = U; du[1] = V; } else // sdim == 3 { const double Jloc[3] = {j(q,0,e), j(q,1,e), j(q,2,e)}; - const double t = 1.0/(Jloc[0]*Jloc[0] + Jloc[1]*Jloc[1] + Jloc[2]*Jloc[2]); - const double U = t*Jloc[0]*du[0]; - const double V = t*Jloc[1]*du[0]; - const double W = t*Jloc[2]*du[0]; + double Jinv[3]; + kernels::CalcLeftInverse<3,1>(Jloc, Jinv); + const double U = Jinv[0]*du[0]; + const double V = Jinv[1]*du[0]; + const double W = Jinv[2]*du[0]; du[0] = U; du[1] = V; du[2] = W; @@ -194,25 +196,10 @@ static void Derivatives2D(const int NE, Jloc[3] = j(qx,qy,0,1,e); Jloc[4] = j(qx,qy,1,1,e); Jloc[5] = j(qx,qy,2,1,e); - - double ee, gg, ff; - ee = Jloc[0]*Jloc[0] + Jloc[1]*Jloc[1] + Jloc[2]*Jloc[2]; - gg = Jloc[3]*Jloc[3] + Jloc[4]*Jloc[4] + Jloc[5]*Jloc[5]; - ff = Jloc[0]*Jloc[3] + Jloc[1]*Jloc[4] + Jloc[2]*Jloc[5]; - const double t = 1.0 / (ee*gg - ff*ff); - ee *= t; gg *= t; ff *= t; - - Jinv[0] = Jloc[0]*gg - Jloc[3]*ff; - Jinv[1] = Jloc[3]*ee - Jloc[0]*ff; - Jinv[2] = Jloc[1]*gg - Jloc[4]*ff; - Jinv[3] = Jloc[4]*ee - Jloc[1]*ff; - Jinv[4] = Jloc[2]*gg - Jloc[5]*ff; - Jinv[5] = Jloc[5]*ee - Jloc[2]*ff; - + kernels::CalcLeftInverse<3,2>(Jloc, Jinv); const double U = Jinv[0]*du[0] + Jinv[1]*du[1]; const double V = Jinv[2]*du[0] + Jinv[3]*du[1]; const double W = Jinv[4]*du[0] + Jinv[5]*du[1]; - du[0] = U; du[1] = V; du[2] = W; diff --git a/linalg/densemat.cpp b/linalg/densemat.cpp index 15ef7c0cfc1..86391f2397f 100644 --- a/linalg/densemat.cpp +++ b/linalg/densemat.cpp @@ -2589,49 +2589,30 @@ void CalcInverse(const DenseMatrix &a, DenseMatrix &inva) MFEM_ASSERT(inva.Height() == a.Width(), "incorrect dimensions"); MFEM_ASSERT(inva.Width() == a.Height(), "incorrect dimensions"); - double t; - if (a.Width() < a.Height()) { const double *d = a.Data(); double *id = inva.Data(); if (a.Height() == 2) { - t = 1.0 / (d[0]*d[0] + d[1]*d[1]); - id[0] = d[0] * t; - id[1] = d[1] * t; + kernels::CalcLeftInverse<2,1>(d, id); } else { if (a.Width() == 1) { - t = 1.0 / (d[0]*d[0] + d[1]*d[1] + d[2]*d[2]); - id[0] = d[0] * t; - id[1] = d[1] * t; - id[2] = d[2] * t; + kernels::CalcLeftInverse<3,1>(d, id); } else { - double e, g, f; - e = d[0]*d[0] + d[1]*d[1] + d[2]*d[2]; - g = d[3]*d[3] + d[4]*d[4] + d[5]*d[5]; - f = d[0]*d[3] + d[1]*d[4] + d[2]*d[5]; - t = 1.0 / (e*g - f*f); - e *= t; g *= t; f *= t; - - id[0] = d[0]*g - d[3]*f; - id[1] = d[3]*e - d[0]*f; - id[2] = d[1]*g - d[4]*f; - id[3] = d[4]*e - d[1]*f; - id[4] = d[2]*g - d[5]*f; - id[5] = d[5]*e - d[2]*f; + kernels::CalcLeftInverse<3,2>(d, id); } } return; } #ifdef MFEM_DEBUG - t = a.Det(); + const double t = a.Det(); MFEM_ASSERT(std::abs(t) > 1.0e-14 * pow(a.FNorm()/a.Width(), a.Width()), "singular matrix!"); #endif diff --git a/linalg/kernels.hpp b/linalg/kernels.hpp index 1c05e84b16f..62230548e02 100644 --- a/linalg/kernels.hpp +++ b/linalg/kernels.hpp @@ -406,6 +406,10 @@ void MultAtB(const int Aheight, const int Awidth, const int Bwidth, } } +/// Given a matrix of size 2x1, 3x1, or 3x2, compute the left inverse. +template MFEM_HOST_DEVICE +void CalcLeftInverse(const double *data, double *left_inv); + /// Compute the spectrum of the matrix of size dim with given @a data, returning /// the eigenvalues in the array @a lambda and the eigenvectors in the array @a /// vec (listed consecutively). @@ -1062,6 +1066,41 @@ int Reduce3S(const int &mode, } // namespace kernels::internal +// Implementations of CalcLeftInverse for dim = 1, 2. + +template<> MFEM_HOST_DEVICE inline +void CalcLeftInverse<2,1>(const double *d, double *left_inv) +{ + const double t = 1.0 / (d[0]*d[0] + d[1]*d[1]); + left_inv[0] = d[0] * t; + left_inv[1] = d[1] * t; +} + +template<> MFEM_HOST_DEVICE inline +void CalcLeftInverse<3,1>(const double *d, double *left_inv) +{ + const double t = 1.0 / (d[0]*d[0] + d[1]*d[1] + d[2]*d[2]); + left_inv[0] = d[0] * t; + left_inv[1] = d[1] * t; + left_inv[2] = d[2] * t; +} + +template<> MFEM_HOST_DEVICE inline +void CalcLeftInverse<3,2>(const double *d, double *left_inv) +{ + double e = d[0]*d[0] + d[1]*d[1] + d[2]*d[2]; + double g = d[3]*d[3] + d[4]*d[4] + d[5]*d[5]; + double f = d[0]*d[3] + d[1]*d[4] + d[2]*d[5]; + const double t = 1.0 / (e*g - f*f); + e *= t; g *= t; f *= t; + + left_inv[0] = d[0]*g - d[3]*f; + left_inv[1] = d[3]*e - d[0]*f; + left_inv[2] = d[1]*g - d[4]*f; + left_inv[3] = d[4]*e - d[1]*f; + left_inv[4] = d[2]*g - d[5]*f; + left_inv[5] = d[5]*e - d[2]*f; +} // Implementations of CalcEigenvalues and CalcSingularvalue for dim = 2, 3.