phynexis-deps/mfem
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge branch 'master' into bugfix/submesh-issue-4107

Browse Source
This commit is contained in:
Tzanio Kolev
2024-02-03 14:27:46 -08:00
committed by GitHub
parent b87027927f 33617ce3a7
commit 1f4ff30101
20 changed files with 121 additions and 71 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fem/bilinearform.cpp
+3 -2
View File
@@ -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
fem/bilininteg.cpp
+5 -5
View File
@@ -2949,16 +2949,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)
fem/coefficient.cpp
+2 -2
View File
@@ -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");
}
fem/coefficient.hpp
+2 -2
View File
@@ -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; }
fem/eltrans.hpp
+7 -4
View File
@@ -170,12 +170,13 @@ 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
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() { }
};
@@ -442,14 +443,16 @@ 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
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);
}
fem/gridfunc.cpp
+9 -8
View File
@@ -2058,7 +2058,7 @@ void GridFunction::AccumulateAndCountZones(VectorCoefficient &vcoeff,
}
void GridFunction::AccumulateAndCountBdrValues(
Coefficient *coeff[], VectorCoefficient *vcoeff, Array<int> &attr,
Coefficient *coeff[], VectorCoefficient *vcoeff, const Array<int> &attr,
Array<int> &values_counter)
{
int i, j, fdof, d, ind, vdim;
@@ -2200,7 +2200,7 @@ static void accumulate_dofs(const Array<int> &dofs, const Vector &vals,
}
void GridFunction::AccumulateAndCountBdrTangentValues(
VectorCoefficient &vcoeff, Array<int> &bdr_attr,
VectorCoefficient &vcoeff, const Array<int> &bdr_attr,
Array<int> &values_counter)
{
const FiniteElement *fe;
@@ -2576,7 +2576,7 @@ void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff,
}
void GridFunction::ProjectBdrCoefficient(VectorCoefficient &vcoeff,
Array<int> &attr)
const Array<int> &attr)
{
Array<int> values_counter;
AccumulateAndCountBdrValues(NULL, &vcoeff, attr, values_counter);
@@ -2593,7 +2593,8 @@ void GridFunction::ProjectBdrCoefficient(VectorCoefficient &vcoeff,
#endif
}
void GridFunction::ProjectBdrCoefficient(Coefficient *coeff[], Array<int> &attr)
void GridFunction::ProjectBdrCoefficient(Coefficient *coeff[],
const Array<int> &attr)
{
Array<int> values_counter;
// this->HostReadWrite(); // done inside the next call
@@ -2623,7 +2624,7 @@ void GridFunction::ProjectBdrCoefficient(Coefficient *coeff[], Array<int> &attr)
}
void GridFunction::ProjectBdrCoefficientNormal(
VectorCoefficient &vcoeff, Array<int> &bdr_attr)
VectorCoefficient &vcoeff, const Array<int> &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<int> &bdr_attr)
VectorCoefficient &vcoeff, const Array<int> &bdr_attr)
{
Array<int> values_counter;
AccumulateAndCountBdrTangentValues(vcoeff, bdr_attr, values_counter);
@@ -2746,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;
@@ -2759,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;
}
fem/gridfunc.hpp
+9 -7
View File
@@ -446,11 +446,12 @@ protected:
Array<int> &zones_per_dof);
void AccumulateAndCountBdrValues(Coefficient *coeff[],
VectorCoefficient *vcoeff, Array<int> &attr,
VectorCoefficient *vcoeff,
const Array<int> &attr,
Array<int> &values_counter);
void AccumulateAndCountBdrTangentValues(VectorCoefficient &vcoeff,
Array<int> &bdr_attr,
const Array<int> &bdr_attr,
Array<int> &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<int> &attr)
void ProjectBdrCoefficient(Coefficient &coeff, const Array<int> &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<int> &attr);
const Array<int> &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<int> &attr);
virtual void ProjectBdrCoefficient(Coefficient *coeff[],
const Array<int> &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<int> &bdr_attr);
const Array<int> &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<int> &bdr_attr);
const Array<int> &bdr_attr);
virtual double ComputeL2Error(Coefficient *exsol[],
const IntegrationRule *irs[] = NULL,
fem/intrules.cpp
+1
View File
@@ -917,6 +917,7 @@ int Quadrature1D::CheckClosed(int type)
{
case GaussLobatto:
case ClosedUniform:
case ClosedGL:
return type;
default:
return Invalid;
fem/lininteg.hpp
+2 -2
View File
@@ -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)
fem/pgridfunc.cpp
+2 -2
View File
@@ -666,7 +666,7 @@ void ParGridFunction::ProjectDiscCoefficient(VectorCoefficient &vcoeff,
}
void ParGridFunction::ProjectBdrCoefficient(
Coefficient *coeff[], VectorCoefficient *vcoeff, Array<int> &attr)
Coefficient *coeff[], VectorCoefficient *vcoeff, const Array<int> &attr)
{
Array<int> values_counter;
AccumulateAndCountBdrValues(coeff, vcoeff, attr, values_counter);
@@ -720,7 +720,7 @@ void ParGridFunction::ProjectBdrCoefficient(
}
void ParGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff,
Array<int> &bdr_attr)
const Array<int> &bdr_attr)
{
Array<int> values_counter;
AccumulateAndCountBdrTangentValues(vcoeff, bdr_attr, values_counter);
fem/pgridfunc.hpp
+5 -4
View File
@@ -44,7 +44,7 @@ protected:
Vector send_data;
void ProjectBdrCoefficient(Coefficient *coeff[], VectorCoefficient *vcoeff,
Array<int> &attr);
const Array<int> &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<int> &attr) override
const Array<int> &attr) override
{ ProjectBdrCoefficient(NULL, &vcoeff, attr); }
// Only the values in the master are guaranteed to be correct!
void ProjectBdrCoefficient(Coefficient *coeff[], Array<int> &attr) override
void ProjectBdrCoefficient(Coefficient *coeff[],
const Array<int> &attr) override
{ ProjectBdrCoefficient(coeff, NULL, attr); }
// Only the values in the master are guaranteed to be correct!
void ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff,
Array<int> &bdr_attr) override;
const Array<int> &bdr_attr) override;
double ComputeL1Error(Coefficient *exsol[],
const IntegrationRule *irs[] = NULL) const override
fem/transfer.cpp
+42 -12
View File
@@ -1420,14 +1420,20 @@ void Prolongation2D(const int NE, const int D1D, const int Q1D,
const Array<double>& 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<double>& 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<double>& 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<double>& 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];
linalg/blockoperator.cpp
+2 -2
View File
@@ -45,12 +45,12 @@ BlockOperator::BlockOperator(const Array<int> & row_offsets_,
op = static_cast<Operator *>(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))
{
linalg/blockoperator.hpp
+3 -6
View File
@@ -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<int> col_offsets;
//! 2D array that stores each block of the operator.
Array2D<Operator *> op;
Array2D<const Operator *> op;
//! 2D array that stores a coefficient for each block of the operator.
Array2D<double> coef;
linalg/operator.cpp
+10 -4
View File
@@ -481,7 +481,8 @@ ConstrainedOperator::ConstrainedOperator(Operator *A, const Array<int> &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<int> &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();
linalg/operator.hpp
+6 -4
View File
@@ -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() override { if (own_A) { delete A; } }
};
/** @brief Rectangular Operator for imposing essential boundary conditions on
linalg/sparsemat.cpp
+1
View File
@@ -1767,6 +1767,7 @@ void SparseMatrix::EliminateCol(int col, DiagonalPolicy dpolicy)
if (aux->Column == col)
{
aux->Value = 0.0;
break;
}
}
}
mesh/pmesh.hpp
+4
View File
@@ -491,6 +491,10 @@ public:
void GenerateOffsets(int N, HYPRE_BigInt loc_sizes[],
Array<HYPRE_BigInt> *offsets[]) const;
/** 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();
miniapps/solvers/div_free_solver.cpp
+5 -4
View File
@@ -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<HypreParMatrix&>(ops_[l]->GetBlock(0, 0));
auto& B_f = static_cast<HypreParMatrix&>(ops_[l]->GetBlock(1, 0));
auto& M_f = static_cast<const HypreParMatrix&>(ops_[l]->GetBlock(0, 0));
auto& B_f = static_cast<const HypreParMatrix&>(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<HypreParMatrix&>(ops_.Last()->GetBlock(0, 0));
auto& M = dynamic_cast<const HypreParMatrix&>(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);
miniapps/solvers/div_free_solver.hpp
+1 -1
View File
@@ -190,7 +190,7 @@ class BDPMinresSolver : public DarcySolver
MINRESSolver solver_;
Array<int> 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) { }