From a8708fc6cae2c826f3908198f97940edbd9998cb Mon Sep 17 00:00:00 2001 From: Will Pazner Date: Wed, 10 Jan 2024 09:43:43 -0800 Subject: [PATCH 01/20] Override AddMult for ConstrainedOperator Use existing workspace vector to avoid extra allocations --- linalg/operator.cpp | 14 ++++++++++---- linalg/operator.hpp | 10 ++++++---- 2 files changed, 16 insertions(+), 8 deletions(-) diff --git a/linalg/operator.cpp b/linalg/operator.cpp index 1f214ece7a..f92f34be6d 100644 --- a/linalg/operator.cpp +++ b/linalg/operator.cpp @@ -481,7 +481,8 @@ ConstrainedOperator::ConstrainedOperator(Operator *A, const Array &list, MemoryType mem_type = GetMemoryType(mem_class); list.Read(); // TODO: just ensure 'list' is registered, no need to copy it constraint_list.MakeRef(list); - // typically z and w are large vectors, so store them on the device + // typically z and w are large vectors, so use the device (GPU) to perform + // operations on them z.SetSize(height, mem_type); z.UseDevice(true); w.SetSize(height, mem_type); w.UseDevice(true); } @@ -591,6 +592,13 @@ void ConstrainedOperator::Mult(const Vector &x, Vector &y) const } } +void ConstrainedOperator::AddMult(const Vector &x, Vector &y, + const double a) const +{ + Mult(x, w); + y.Add(a, w); +} + RectangularConstrainedOperator::RectangularConstrainedOperator( Operator *A, const Array &trial_list, @@ -625,9 +633,7 @@ void RectangularConstrainedOperator::EliminateRHS(const Vector &x, d_w[id] = d_x[id]; }); - // A.AddMult(w, b, -1.0); // if available to all Operators - A->Mult(w, z); - b -= z; + A->AddMult(w, b, -1.0); const int test_csz = test_constraints.Size(); auto test_idx = test_constraints.Read(); diff --git a/linalg/operator.hpp b/linalg/operator.hpp index baa9bf7672..30b9298f77 100644 --- a/linalg/operator.hpp +++ b/linalg/operator.hpp @@ -892,14 +892,14 @@ public: DiagonalPolicy diag_policy = DIAG_ONE); /// Returns the type of memory in which the solution and temporaries are stored. - virtual MemoryClass GetMemoryClass() const { return mem_class; } + MemoryClass GetMemoryClass() const override { return mem_class; } /// Set the diagonal policy for the constrained operator. void SetDiagonalPolicy(const DiagonalPolicy diag_policy_) { diag_policy = diag_policy_; } /// Diagonal of A, modified according to the used DiagonalPolicy. - virtual void AssembleDiagonal(Vector &diag) const; + void AssembleDiagonal(Vector &diag) const override; /** @brief Eliminate "essential boundary condition" values specified in @a x from the given right-hand side @a b. @@ -922,10 +922,12 @@ public: where the "_b" subscripts denote the essential (boundary) indices/dofs of the vectors, and "_i" -- the rest of the entries. */ - virtual void Mult(const Vector &x, Vector &y) const; + void Mult(const Vector &x, Vector &y) const override; + + void AddMult(const Vector &x, Vector &y, const double a = 1.0) const override; /// Destructor: destroys the unconstrained Operator, if owned. - virtual ~ConstrainedOperator() { if (own_A) { delete A; } } + ~ConstrainedOperator() { if (own_A) { delete A; } } }; /** @brief Rectangular Operator for imposing essential boundary conditions on From 68e513c40eed5bab68bbcc6820eecb7f9de092dd Mon Sep 17 00:00:00 2001 From: Will Pazner Date: Thu, 11 Jan 2024 10:06:28 -0800 Subject: [PATCH 02/20] Use Write instead of ReadWrite in TensorProductPRefinementTransferOperator Can avoid some unnecessary H to D memcpy --- fem/transfer.cpp | 54 +++++++++++++++++++++++++++++++++++++----------- 1 file changed, 42 insertions(+), 12 deletions(-) diff --git a/fem/transfer.cpp b/fem/transfer.cpp index 7c08f5efff..7e3ff4a016 100644 --- a/fem/transfer.cpp +++ b/fem/transfer.cpp @@ -1420,14 +1420,20 @@ void Prolongation2D(const int NE, const int D1D, const int Q1D, const Array& B, const Vector& mask) { auto x_ = Reshape(localL.Read(), D1D, D1D, NE); - auto y_ = Reshape(localH.ReadWrite(), Q1D, Q1D, NE); + auto y_ = Reshape(localH.Write(), Q1D, Q1D, NE); auto B_ = Reshape(B.Read(), Q1D, D1D); auto m_ = Reshape(mask.Read(), Q1D, Q1D, NE); - localH = 0.0; - mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e) { + for (int qy = 0; qy < Q1D; ++qy) + { + for (int qx = 0; qx < Q1D; ++qx) + { + y_(qx, qy, e) = 0.0; + } + } + for (int dy = 0; dy < D1D; ++dy) { double sol_x[DofQuadLimits::MAX_Q1D]; @@ -1467,14 +1473,23 @@ void Prolongation3D(const int NE, const int D1D, const int Q1D, const Array& B, const Vector& mask) { auto x_ = Reshape(localL.Read(), D1D, D1D, D1D, NE); - auto y_ = Reshape(localH.ReadWrite(), Q1D, Q1D, Q1D, NE); + auto y_ = Reshape(localH.Write(), Q1D, Q1D, Q1D, NE); auto B_ = Reshape(B.Read(), Q1D, D1D); auto m_ = Reshape(mask.Read(), Q1D, Q1D, Q1D, NE); - localH = 0.0; - mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e) { + for (int qz = 0; qz < Q1D; ++qz) + { + for (int qy = 0; qy < Q1D; ++qy) + { + for (int qx = 0; qx < Q1D; ++qx) + { + y_(qx, qy, qz, e) = 0.0; + } + } + } + for (int dz = 0; dz < D1D; ++dz) { double sol_xy[DofQuadLimits::MAX_Q1D][DofQuadLimits::MAX_Q1D]; @@ -1539,14 +1554,20 @@ void Restriction2D(const int NE, const int D1D, const int Q1D, const Array& Bt, const Vector& mask) { auto x_ = Reshape(localH.Read(), Q1D, Q1D, NE); - auto y_ = Reshape(localL.ReadWrite(), D1D, D1D, NE); + auto y_ = Reshape(localL.Write(), D1D, D1D, NE); auto Bt_ = Reshape(Bt.Read(), D1D, Q1D); auto m_ = Reshape(mask.Read(), Q1D, Q1D, NE); - localL = 0.0; - mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e) { + for (int dy = 0; dy < D1D; ++dy) + { + for (int dx = 0; dx < D1D; ++dx) + { + y_(dx, dy, e) = 0.0; + } + } + for (int qy = 0; qy < Q1D; ++qy) { double sol_x[DofQuadLimits::MAX_D1D]; @@ -1578,14 +1599,23 @@ void Restriction3D(const int NE, const int D1D, const int Q1D, const Array& Bt, const Vector& mask) { auto x_ = Reshape(localH.Read(), Q1D, Q1D, Q1D, NE); - auto y_ = Reshape(localL.ReadWrite(), D1D, D1D, D1D, NE); + auto y_ = Reshape(localL.Write(), D1D, D1D, D1D, NE); auto Bt_ = Reshape(Bt.Read(), D1D, Q1D); auto m_ = Reshape(mask.Read(), Q1D, Q1D, Q1D, NE); - localL = 0.0; - mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e) { + for (int dz = 0; dz < D1D; ++dz) + { + for (int dy = 0; dy < D1D; ++dy) + { + for (int dx = 0; dx < D1D; ++dx) + { + y_(dx, dy, dz, e) = 0.0; + } + } + } + for (int qz = 0; qz < Q1D; ++qz) { double sol_xy[DofQuadLimits::MAX_D1D][DofQuadLimits::MAX_D1D]; From c54bc0aa55c6b9823a98b9f294ba9bd9323ce763 Mon Sep 17 00:00:00 2001 From: Robert Carson Date: Thu, 11 Jan 2024 10:48:14 -0800 Subject: [PATCH 03/20] Add ability to set tolerance in TransformBack func of the ElementTransformation class --- fem/eltrans.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/fem/eltrans.hpp b/fem/eltrans.hpp index 8d3955ab41..57d259cb42 100644 --- a/fem/eltrans.hpp +++ b/fem/eltrans.hpp @@ -175,7 +175,7 @@ public: point in physical space. If the inversion fails a non-zero value is returned. This method is not 100 percent reliable for non-linear transformations. */ - virtual int TransformBack(const Vector &pt, IntegrationPoint &ip) = 0; + virtual int TransformBack(const Vector &pt, IntegrationPoint &ip, const double phys_tol = 1e-15) = 0; virtual ~ElementTransformation() { } }; @@ -447,9 +447,10 @@ public: point in physical space. If the inversion fails a non-zero value is returned. This method is not 100 percent reliable for non-linear transformations. */ - virtual int TransformBack(const Vector & v, IntegrationPoint & ip) + virtual int TransformBack(const Vector & v, IntegrationPoint & ip, const double phys_rel_tol = 1e-15) { InverseElementTransformation inv_tr(this); + inv_tr.SetPhysicalRelTol(phys_rel_tol); return inv_tr.Transform(v, ip); } From 6c4f179751ecfa45863a376a805878882d90afb7 Mon Sep 17 00:00:00 2001 From: Robert Carson Date: Thu, 11 Jan 2024 10:56:53 -0800 Subject: [PATCH 04/20] make style --- fem/eltrans.hpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/fem/eltrans.hpp b/fem/eltrans.hpp index 57d259cb42..a81952ad95 100644 --- a/fem/eltrans.hpp +++ b/fem/eltrans.hpp @@ -175,7 +175,8 @@ public: point in physical space. If the inversion fails a non-zero value is returned. This method is not 100 percent reliable for non-linear transformations. */ - virtual int TransformBack(const Vector &pt, IntegrationPoint &ip, const double phys_tol = 1e-15) = 0; + virtual int TransformBack(const Vector &pt, IntegrationPoint &ip, + const double phys_tol = 1e-15) = 0; virtual ~ElementTransformation() { } }; @@ -447,7 +448,8 @@ public: point in physical space. If the inversion fails a non-zero value is returned. This method is not 100 percent reliable for non-linear transformations. */ - virtual int TransformBack(const Vector & v, IntegrationPoint & ip, const double phys_rel_tol = 1e-15) + virtual int TransformBack(const Vector & v, IntegrationPoint & ip, + const double phys_rel_tol = 1e-15) { InverseElementTransformation inv_tr(this); inv_tr.SetPhysicalRelTol(phys_rel_tol); From 4c547f52c247a39618b1be68e58832ac87578ca2 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Mon, 23 Dec 2019 10:28:16 +0100 Subject: [PATCH 05/20] Minor optimization of SparseMatrix::EliminateCol(). --- linalg/sparsemat.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/linalg/sparsemat.cpp b/linalg/sparsemat.cpp index 6f1b0d6828..9ac3b79951 100644 --- a/linalg/sparsemat.cpp +++ b/linalg/sparsemat.cpp @@ -1767,6 +1767,7 @@ void SparseMatrix::EliminateCol(int col, DiagonalPolicy dpolicy) if (aux->Column == col) { aux->Value = 0.0; + break; } } } From d843cd6e8a65f6266383bbf4a317b08b0dad40d7 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Sun, 29 Dec 2019 21:29:18 +0100 Subject: [PATCH 06/20] Added const qualifiers to Operator arguments of BlockOperator. --- linalg/blockoperator.cpp | 4 ++-- linalg/blockoperator.hpp | 9 +++------ 2 files changed, 5 insertions(+), 8 deletions(-) diff --git a/linalg/blockoperator.cpp b/linalg/blockoperator.cpp index b1c0c3c311..c53d8dc36c 100644 --- a/linalg/blockoperator.cpp +++ b/linalg/blockoperator.cpp @@ -45,12 +45,12 @@ BlockOperator::BlockOperator(const Array & row_offsets_, op = static_cast(NULL); } -void BlockOperator::SetDiagonalBlock(int iblock, Operator *opt, double c) +void BlockOperator::SetDiagonalBlock(int iblock, const Operator *opt, double c) { SetBlock(iblock, iblock, opt, c); } -void BlockOperator::SetBlock(int iRow, int iCol, Operator *opt, double c) +void BlockOperator::SetBlock(int iRow, int iCol, const Operator *opt, double c) { if (owns_blocks && op(iRow, iCol)) { diff --git a/linalg/blockoperator.hpp b/linalg/blockoperator.hpp index c9f9c8e46a..8c22cfe7a9 100644 --- a/linalg/blockoperator.hpp +++ b/linalg/blockoperator.hpp @@ -61,14 +61,14 @@ public: * op: the Operator to be inserted. * c: optional scalar multiple for this block. */ - void SetDiagonalBlock(int iblock, Operator *op, double c = 1.0); + void SetDiagonalBlock(int iblock, const Operator *op, double c = 1.0); //! Add a block op in the block-entry (iblock, jblock). /** * irow, icol: The block will be inserted in location (irow, icol). * op: the Operator to be inserted. * c: optional scalar multiple for this block. */ - void SetBlock(int iRow, int iCol, Operator *op, double c = 1.0); + void SetBlock(int iRow, int iCol, const Operator *op, double c = 1.0); //! Return the number of row blocks int NumRowBlocks() const { return nRowBlocks; } @@ -78,9 +78,6 @@ public: //! Check if block (i,j) is a zero block int IsZeroBlock(int i, int j) const { return (op(i,j)==NULL) ? 1 : 0; } //! Return a reference to block i,j - Operator & GetBlock(int i, int j) - { MFEM_VERIFY(op(i,j), ""); return *op(i,j); } - //! Return a reference to block i,j (const version) const Operator & GetBlock(int i, int j) const { MFEM_VERIFY(op(i,j), ""); return *op(i,j); } //! Return the coefficient for block i,j @@ -123,7 +120,7 @@ private: //! Column offsets for the starting position of each block Array col_offsets; //! 2D array that stores each block of the operator. - Array2D op; + Array2D op; //! 2D array that stores a coefficient for each block of the operator. Array2D coef; From f8040db3ec84ac2017a02294a65087e62520197e Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Fri, 24 Dec 2021 08:36:53 +0100 Subject: [PATCH 07/20] Made available FaceIsTrueInterior in ParMesh as it is supposed to be probably. --- mesh/pmesh.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/mesh/pmesh.hpp b/mesh/pmesh.hpp index ce5dbfbc28..48a74a763a 100644 --- a/mesh/pmesh.hpp +++ b/mesh/pmesh.hpp @@ -491,6 +491,8 @@ public: void GenerateOffsets(int N, HYPRE_BigInt loc_sizes[], Array *offsets[]) const; + using Mesh::FaceIsTrueInterior; + void ExchangeFaceNbrData(); void ExchangeFaceNbrNodes(); From 755eb9d65f9ec8061453ce82c462aef435698496 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Wed, 8 Jun 2022 09:07:30 +0200 Subject: [PATCH 08/20] Added const qualifier to quadrature function coefficient. --- fem/coefficient.cpp | 4 ++-- fem/coefficient.hpp | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/fem/coefficient.cpp b/fem/coefficient.cpp index 90eef2c9ce..cebdf348f5 100644 --- a/fem/coefficient.cpp +++ b/fem/coefficient.cpp @@ -1568,7 +1568,7 @@ double ComputeGlobalLpNorm(double p, VectorCoefficient &coeff, ParMesh &pmesh, #endif VectorQuadratureFunctionCoefficient::VectorQuadratureFunctionCoefficient( - QuadratureFunction &qf) + const QuadratureFunction &qf) : VectorCoefficient(qf.GetVDim()), QuadF(qf), index(0) { } void VectorQuadratureFunctionCoefficient::SetComponent(int index_, int length_) @@ -1622,7 +1622,7 @@ void VectorQuadratureFunctionCoefficient::Project(QuadratureFunction &qf) } QuadratureFunctionCoefficient::QuadratureFunctionCoefficient( - QuadratureFunction &qf) : QuadF(qf) + const QuadratureFunction &qf) : QuadF(qf) { MFEM_VERIFY(qf.GetVDim() == 1, "QuadratureFunction's vdim must be 1"); } diff --git a/fem/coefficient.hpp b/fem/coefficient.hpp index 568c9376f0..a80f6e81c6 100644 --- a/fem/coefficient.hpp +++ b/fem/coefficient.hpp @@ -2245,7 +2245,7 @@ private: public: /// Constructor with a quadrature function as input - VectorQuadratureFunctionCoefficient(QuadratureFunction &qf); + VectorQuadratureFunctionCoefficient(const QuadratureFunction &qf); /** Set the starting index within the QuadFunc that'll be used to project outwards as well as the corresponding length. The projected length should @@ -2273,7 +2273,7 @@ private: public: /// Constructor with a quadrature function as input - QuadratureFunctionCoefficient(QuadratureFunction &qf); + QuadratureFunctionCoefficient(const QuadratureFunction &qf); const QuadratureFunction& GetQuadFunction() const { return QuadF; } From 34aac72833f81decc783221a4848cbc2de5f3775 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Wed, 8 Jun 2022 09:08:47 +0200 Subject: [PATCH 09/20] Added default NULL value of the integration rule for quadrature integrators as it is not used anyway. --- fem/lininteg.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/fem/lininteg.hpp b/fem/lininteg.hpp index 0cc5a80d44..25741e0626 100644 --- a/fem/lininteg.hpp +++ b/fem/lininteg.hpp @@ -728,7 +728,7 @@ private: public: VectorQuadratureLFIntegrator(VectorQuadratureFunctionCoefficient &vqfc, - const IntegrationRule *ir) + const IntegrationRule *ir = NULL) : LinearFormIntegrator(ir), vqfc(vqfc) { if (ir) @@ -760,7 +760,7 @@ private: public: QuadratureLFIntegrator(QuadratureFunctionCoefficient &qfc, - const IntegrationRule *ir) + const IntegrationRule *ir = NULL) : LinearFormIntegrator(ir), qfc(qfc) { if (ir) From e311b0034937cde09de42e4a23bbe746661ce21f Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Wed, 13 Jul 2022 11:58:12 +0200 Subject: [PATCH 10/20] Added closed GL quadrature into the closed ones. --- fem/intrules.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/fem/intrules.cpp b/fem/intrules.cpp index cb9544852f..ddce39e0a4 100644 --- a/fem/intrules.cpp +++ b/fem/intrules.cpp @@ -917,6 +917,7 @@ int Quadrature1D::CheckClosed(int type) { case GaussLobatto: case ClosedUniform: + case ClosedGL: return type; default: return Invalid; From 55a35c881814c3708e3a34acdc44b417b79b8043 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Fri, 10 Feb 2023 10:04:10 +0100 Subject: [PATCH 11/20] Fixed boundary trace integration in MixedBilinearForm. --- fem/bilinearform.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/fem/bilinearform.cpp b/fem/bilinearform.cpp index 2661ff7f13..5d08386690 100644 --- a/fem/bilinearform.cpp +++ b/fem/bilinearform.cpp @@ -1585,9 +1585,10 @@ void MixedBilinearForm::Assemble(int skip_zeros) ftr = mesh->GetBdrFaceTransformations(i); if (ftr) { - trial_fes->GetFaceVDofs(ftr->ElementNo, trial_vdofs); + const int iface = mesh->GetBdrElementFaceIndex(i); + trial_fes->GetFaceVDofs(iface, trial_vdofs); test_fes->GetElementVDofs(ftr->Elem1No, test_vdofs); - trial_face_fe = trial_fes->GetFaceElement(ftr->ElementNo); + trial_face_fe = trial_fes->GetFaceElement(iface); test_fe1 = test_fes->GetFE(ftr->Elem1No); // The test_fe2 object is really a dummy and not used on the // boundaries, but we can't dereference a NULL pointer, and we don't From 63c97160564286cab93902b93aebc8c877f14117 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Wed, 15 Feb 2023 16:52:15 +0100 Subject: [PATCH 12/20] Fixed VectorDiffusionOperator::AssembleElementVector(). --- fem/bilininteg.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/fem/bilininteg.cpp b/fem/bilininteg.cpp index 59a35f4606..f7f253851d 100644 --- a/fem/bilininteg.cpp +++ b/fem/bilininteg.cpp @@ -2948,16 +2948,16 @@ void VectorDiffusionIntegrator::AssembleElementVector( } dshape.SetSize(dof, dim); - dshapedxt.SetSize(dof, dim); - // pelmat.SetSize(dim); + dshapedxt.SetSize(dof, sdim); + pelmat.SetSize(dof); - elvect.SetSize(dim*dof); + elvect.SetSize(vdim*dof); // NOTE: DenseMatrix is in column-major order. This is consistent with // vectors ordered byNODES. In the resulting DenseMatrix, each column // corresponds to a particular vdim. - DenseMatrix mat_in(elfun.GetData(), dof, dim); - DenseMatrix mat_out(elvect.GetData(), dof, dim); + DenseMatrix mat_in(elfun.GetData(), dof, vdim); + DenseMatrix mat_out(elvect.GetData(), dof, vdim); const IntegrationRule *ir = IntRule; if (ir == NULL) From 2d582d80c9c7e605c1cc9afe8a50647687bd8f84 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Tue, 7 Mar 2023 09:32:27 +0100 Subject: [PATCH 13/20] Added const qualifier for boundary attributes in GridFunction projection methods. --- fem/gridfunc.cpp | 13 +++++++------ fem/gridfunc.hpp | 16 +++++++++------- 2 files changed, 16 insertions(+), 13 deletions(-) diff --git a/fem/gridfunc.cpp b/fem/gridfunc.cpp index 310d8d7043..6f2a930586 100644 --- a/fem/gridfunc.cpp +++ b/fem/gridfunc.cpp @@ -2058,7 +2058,7 @@ void GridFunction::AccumulateAndCountZones(VectorCoefficient &vcoeff, } void GridFunction::AccumulateAndCountBdrValues( - Coefficient *coeff[], VectorCoefficient *vcoeff, Array &attr, + Coefficient *coeff[], VectorCoefficient *vcoeff, const Array &attr, Array &values_counter) { int i, j, fdof, d, ind, vdim; @@ -2200,7 +2200,7 @@ static void accumulate_dofs(const Array &dofs, const Vector &vals, } void GridFunction::AccumulateAndCountBdrTangentValues( - VectorCoefficient &vcoeff, Array &bdr_attr, + VectorCoefficient &vcoeff, const Array &bdr_attr, Array &values_counter) { const FiniteElement *fe; @@ -2576,7 +2576,7 @@ void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff, } void GridFunction::ProjectBdrCoefficient(VectorCoefficient &vcoeff, - Array &attr) + const Array &attr) { Array values_counter; AccumulateAndCountBdrValues(NULL, &vcoeff, attr, values_counter); @@ -2593,7 +2593,8 @@ void GridFunction::ProjectBdrCoefficient(VectorCoefficient &vcoeff, #endif } -void GridFunction::ProjectBdrCoefficient(Coefficient *coeff[], Array &attr) +void GridFunction::ProjectBdrCoefficient(Coefficient *coeff[], + const Array &attr) { Array values_counter; // this->HostReadWrite(); // done inside the next call @@ -2623,7 +2624,7 @@ void GridFunction::ProjectBdrCoefficient(Coefficient *coeff[], Array &attr) } void GridFunction::ProjectBdrCoefficientNormal( - VectorCoefficient &vcoeff, Array &bdr_attr) + VectorCoefficient &vcoeff, const Array &bdr_attr) { #if 0 // implementation for the case when the face dofs are integrals of the @@ -2698,7 +2699,7 @@ void GridFunction::ProjectBdrCoefficientNormal( } void GridFunction::ProjectBdrCoefficientTangent( - VectorCoefficient &vcoeff, Array &bdr_attr) + VectorCoefficient &vcoeff, const Array &bdr_attr) { Array values_counter; AccumulateAndCountBdrTangentValues(vcoeff, bdr_attr, values_counter); diff --git a/fem/gridfunc.hpp b/fem/gridfunc.hpp index d7cf303be5..e6330c6a1c 100644 --- a/fem/gridfunc.hpp +++ b/fem/gridfunc.hpp @@ -446,11 +446,12 @@ protected: Array &zones_per_dof); void AccumulateAndCountBdrValues(Coefficient *coeff[], - VectorCoefficient *vcoeff, Array &attr, + VectorCoefficient *vcoeff, + const Array &attr, Array &values_counter); void AccumulateAndCountBdrTangentValues(VectorCoefficient &vcoeff, - Array &bdr_attr, + const Array &bdr_attr, Array &values_counter); // Complete the computation of averages; called e.g. after @@ -465,7 +466,7 @@ public: /** @brief Project a Coefficient on the GridFunction, modifying only DOFs on the boundary associated with the boundary attributes marked in the @a attr array. */ - void ProjectBdrCoefficient(Coefficient &coeff, Array &attr) + void ProjectBdrCoefficient(Coefficient &coeff, const Array &attr) { Coefficient *coeff_p = &coeff; ProjectBdrCoefficient(&coeff_p, attr); @@ -475,26 +476,27 @@ public: DOFs on the boundary associated with the boundary attributes marked in the @a attr array. */ virtual void ProjectBdrCoefficient(VectorCoefficient &vcoeff, - Array &attr); + const Array &attr); /** @brief Project a set of Coefficient%s on the components of the GridFunction, modifying only DOFs on the boundary associated with the boundary attributed marked in the @a attr array. */ /** If a Coefficient pointer in the array @a coeff is NULL, that component will not be touched. */ - virtual void ProjectBdrCoefficient(Coefficient *coeff[], Array &attr); + virtual void ProjectBdrCoefficient(Coefficient *coeff[], + const Array &attr); /** Project the normal component of the given VectorCoefficient on the boundary. Only boundary attributes that are marked in 'bdr_attr' are projected. Assumes RT-type VectorFE GridFunction. */ void ProjectBdrCoefficientNormal(VectorCoefficient &vcoeff, - Array &bdr_attr); + const Array &bdr_attr); /** @brief Project the tangential components of the given VectorCoefficient on the boundary. Only boundary attributes that are marked in @a bdr_attr are projected. Assumes ND-type VectorFE GridFunction. */ virtual void ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff, - Array &bdr_attr); + const Array &bdr_attr); virtual double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[] = NULL, From 20dbb8efbae48fe23c78887d46d23d22b1a45cf1 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 7 Sep 2023 10:33:58 +0200 Subject: [PATCH 14/20] Added support of integral elements to L2 error. --- fem/gridfunc.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/fem/gridfunc.cpp b/fem/gridfunc.cpp index 6f2a930586..0d719935fb 100644 --- a/fem/gridfunc.cpp +++ b/fem/gridfunc.cpp @@ -2747,7 +2747,8 @@ double GridFunction::ComputeL2Error( for (j = 0; j < ir->GetNPoints(); j++) { const IntegrationPoint &ip = ir->IntPoint(j); - fe->CalcShape(ip, shape); + transf->SetIntPoint(&ip); + fe->CalcPhysShape(*transf, shape); for (d = 0; d < fes->GetVDim(); d++) { a = 0; @@ -2760,7 +2761,6 @@ double GridFunction::ComputeL2Error( { a -= (*this)(-1-vdofs[fdof*d+k]) * shape(k); } - transf->SetIntPoint(&ip); a -= exsol[d]->Eval(*transf, ip); error += ip.weight * transf->Weight() * a * a; } From a0ae98f4001204ccf9f4a731c09b95191f8441b7 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Mon, 22 Jan 2024 16:23:43 -0800 Subject: [PATCH 15/20] Added const qualifier to boundary attributes in ParGridFunction projection methods. --- fem/pgridfunc.cpp | 4 ++-- fem/pgridfunc.hpp | 9 +++++---- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/fem/pgridfunc.cpp b/fem/pgridfunc.cpp index e546093599..178d5fc5cc 100644 --- a/fem/pgridfunc.cpp +++ b/fem/pgridfunc.cpp @@ -666,7 +666,7 @@ void ParGridFunction::ProjectDiscCoefficient(VectorCoefficient &vcoeff, } void ParGridFunction::ProjectBdrCoefficient( - Coefficient *coeff[], VectorCoefficient *vcoeff, Array &attr) + Coefficient *coeff[], VectorCoefficient *vcoeff, const Array &attr) { Array values_counter; AccumulateAndCountBdrValues(coeff, vcoeff, attr, values_counter); @@ -720,7 +720,7 @@ void ParGridFunction::ProjectBdrCoefficient( } void ParGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff, - Array &bdr_attr) + const Array &bdr_attr) { Array values_counter; AccumulateAndCountBdrTangentValues(vcoeff, bdr_attr, values_counter); diff --git a/fem/pgridfunc.hpp b/fem/pgridfunc.hpp index 58c6c02862..6f2ec89246 100644 --- a/fem/pgridfunc.hpp +++ b/fem/pgridfunc.hpp @@ -44,7 +44,7 @@ protected: Vector send_data; void ProjectBdrCoefficient(Coefficient *coeff[], VectorCoefficient *vcoeff, - Array &attr); + const Array &attr); public: ParGridFunction() { pfes = NULL; } @@ -256,16 +256,17 @@ public: // Only the values in the master are guaranteed to be correct! void ProjectBdrCoefficient(VectorCoefficient &vcoeff, - Array &attr) override + const Array &attr) override { ProjectBdrCoefficient(NULL, &vcoeff, attr); } // Only the values in the master are guaranteed to be correct! - void ProjectBdrCoefficient(Coefficient *coeff[], Array &attr) override + void ProjectBdrCoefficient(Coefficient *coeff[], + const Array &attr) override { ProjectBdrCoefficient(coeff, NULL, attr); } // Only the values in the master are guaranteed to be correct! void ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff, - Array &bdr_attr) override; + const Array &bdr_attr) override; double ComputeL1Error(Coefficient *exsol[], const IntegrationRule *irs[] = NULL) const override From 7a2ab45ae02c2bfd7b02edd99ce28055cb743354 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Mon, 22 Jan 2024 16:55:41 -0800 Subject: [PATCH 16/20] Fixed div-free solver miniapp to work with constant operators. --- miniapps/solvers/div_free_solver.cpp | 9 +++++---- miniapps/solvers/div_free_solver.hpp | 2 +- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/miniapps/solvers/div_free_solver.cpp b/miniapps/solvers/div_free_solver.cpp index 978abbcc40..7c3905e061 100644 --- a/miniapps/solvers/div_free_solver.cpp +++ b/miniapps/solvers/div_free_solver.cpp @@ -314,7 +314,8 @@ void SaddleSchwarzSmoother::Mult(const Vector & x, Vector & y) const blk_y.GetBlock(1) -= coarse_l2_projection; } -BDPMinresSolver::BDPMinresSolver(HypreParMatrix& M, HypreParMatrix& B, +BDPMinresSolver::BDPMinresSolver(const HypreParMatrix& M, + const HypreParMatrix& B, IterSolveParameters param) : DarcySolver(M.NumRows(), B.NumRows()), op_(offsets_), prec_(offsets_), BT_(B.Transpose()), solver_(M.GetComm()) @@ -360,8 +361,8 @@ DivFreeSolver::DivFreeSolver(const HypreParMatrix &M, const HypreParMatrix& B, for (int l = data.P_l2.size(); l >= 0; --l) { - auto& M_f = static_cast(ops_[l]->GetBlock(0, 0)); - auto& B_f = static_cast(ops_[l]->GetBlock(1, 0)); + auto& M_f = static_cast(ops_[l]->GetBlock(0, 0)); + auto& B_f = static_cast(ops_[l]->GetBlock(1, 0)); if (l == 0) { @@ -586,7 +587,7 @@ void DivFreeSolver::Mult(const Vector & x, Vector & y) const ch.Clear(); ch.Start(); - auto M = dynamic_cast(ops_.Last()->GetBlock(0, 0)); + auto M = dynamic_cast(ops_.Last()->GetBlock(0, 0)); M.Mult(-1.0, correction.GetBlock(0), 1.0, resid.GetBlock(0)); SolvePotential(resid.GetBlock(0), correction.GetBlock(1)); blk_y.GetBlock(1) += correction.GetBlock(1); diff --git a/miniapps/solvers/div_free_solver.hpp b/miniapps/solvers/div_free_solver.hpp index ad16c2c45a..57177231fc 100644 --- a/miniapps/solvers/div_free_solver.hpp +++ b/miniapps/solvers/div_free_solver.hpp @@ -190,7 +190,7 @@ class BDPMinresSolver : public DarcySolver MINRESSolver solver_; Array ess_zero_dofs_; public: - BDPMinresSolver(HypreParMatrix& M, HypreParMatrix& B, + BDPMinresSolver(const HypreParMatrix& M, const HypreParMatrix& B, IterSolveParameters param); virtual void Mult(const Vector & x, Vector & y) const; virtual void SetOperator(const Operator &op) { } From c3a455d86f06f2e2cde7a536da6de4c8183a1658 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Wed, 24 Jan 2024 12:23:20 -0800 Subject: [PATCH 17/20] Fixed auto reference in div_free_solver. --- miniapps/solvers/div_free_solver.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/miniapps/solvers/div_free_solver.cpp b/miniapps/solvers/div_free_solver.cpp index 7c3905e061..e646b3bc20 100644 --- a/miniapps/solvers/div_free_solver.cpp +++ b/miniapps/solvers/div_free_solver.cpp @@ -587,7 +587,7 @@ void DivFreeSolver::Mult(const Vector & x, Vector & y) const ch.Clear(); ch.Start(); - auto M = dynamic_cast(ops_.Last()->GetBlock(0, 0)); + auto& M = dynamic_cast(ops_.Last()->GetBlock(0, 0)); M.Mult(-1.0, correction.GetBlock(0), 1.0, resid.GetBlock(0)); SolvePotential(resid.GetBlock(0), correction.GetBlock(1)); blk_y.GetBlock(1) += correction.GetBlock(1); From 3ae6ec61e5edfddc8bf72c953c50909415f42812 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 25 Jan 2024 08:59:41 -0800 Subject: [PATCH 18/20] Changed using to stub for ParMesh::FaceIsTrueInterior(). --- mesh/pmesh.hpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/mesh/pmesh.hpp b/mesh/pmesh.hpp index 48a74a763a..ddcfd0934a 100644 --- a/mesh/pmesh.hpp +++ b/mesh/pmesh.hpp @@ -491,7 +491,9 @@ public: void GenerateOffsets(int N, HYPRE_BigInt loc_sizes[], Array *offsets[]) const; - using Mesh::FaceIsTrueInterior; + /** Return true if the face is interior or shared. In parallel, this + method only works if the face neighbor data is exchanged. */ + inline bool FaceIsTrueInterior(int FaceNo) const { return Mesh::FaceIsTrueInterior(FaceNo); } void ExchangeFaceNbrData(); void ExchangeFaceNbrNodes(); From 1deac8c6bffef86ec4573e1af96a4da9bc77e2b5 Mon Sep 17 00:00:00 2001 From: Robert Carson Date: Tue, 30 Jan 2024 17:30:28 -0800 Subject: [PATCH 19/20] Add documentation for new param in ElementTransformation::TransformBack --- fem/eltrans.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/fem/eltrans.hpp b/fem/eltrans.hpp index a81952ad95..ca3892d4c2 100644 --- a/fem/eltrans.hpp +++ b/fem/eltrans.hpp @@ -170,7 +170,7 @@ public: virtual int GetSpaceDim() const = 0; /** @brief Transform a point @a pt from physical space to a point @a ip in - reference space. */ + reference space and optionally can set a solver tolerance using @a phys_tol. */ /** Attempt to find the IntegrationPoint that is transformed into the given point in physical space. If the inversion fails a non-zero value is returned. This method is not 100 percent reliable for non-linear @@ -443,7 +443,7 @@ public: virtual int GetSpaceDim() const { return PointMat.Height(); } /** @brief Transform a point @a pt from physical space to a point @a ip in - reference space. */ + reference space and optionally can set a solver tolerance using @a phys_tol. */ /** Attempt to find the IntegrationPoint that is transformed into the given point in physical space. If the inversion fails a non-zero value is returned. This method is not 100 percent reliable for non-linear From 58319e5f94a9a33994aff29ef793d5ea5e6765a2 Mon Sep 17 00:00:00 2001 From: Will Pazner Date: Thu, 1 Feb 2024 06:48:41 -0800 Subject: [PATCH 20/20] Mark ConstrainedOperator::~ConstrainedOperator as override --- linalg/operator.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/linalg/operator.hpp b/linalg/operator.hpp index 3e20b46fb3..aeb3564869 100644 --- a/linalg/operator.hpp +++ b/linalg/operator.hpp @@ -927,7 +927,7 @@ public: void AddMult(const Vector &x, Vector &y, const double a = 1.0) const override; /// Destructor: destroys the unconstrained Operator, if owned. - ~ConstrainedOperator() { if (own_A) { delete A; } } + ~ConstrainedOperator() override { if (own_A) { delete A; } } }; /** @brief Rectangular Operator for imposing essential boundary conditions on