26 lagrange_tdofs(lagrange_tdofs_),
27 primary_tdofs(primary_tdofs_),
28 secondary_tdofs(secondary_tdofs_)
30 MFEM_VERIFY(lagrange_tdofs.
Size() == secondary_tdofs.
Size(),
31 "Dof sizes don't match!");
92 eliminators(eliminators_)
98 MFEM_ASSERT(vin.
Size() ==
width,
"Wrong vector size!");
99 MFEM_ASSERT(vout.
Size() ==
height,
"Wrong vector size!");
103 for (
int k = 0; k < eliminators.Size(); ++k)
116 MFEM_ASSERT(vin.
Size() ==
height,
"Wrong vector size!");
117 MFEM_ASSERT(vout.
Size() ==
width,
"Wrong vector size!");
121 for (
int k = 0; k < eliminators.Size(); ++k)
137 for (
int i = 0; i <
height; ++i)
142 for (
int k = 0; k < eliminators.Size(); ++k)
150 for (
int jz = 0; jz < elim->
PrimaryDofs().Size(); ++jz)
153 mat->
Add(i, j, mat_k(iz, jz));
165 MFEM_ASSERT(rtilde.
Size() == Aop.
Height(),
"Sizes don't match!");
168 for (
int k = 0; k < eliminators.Size(); ++k)
183 MFEM_ASSERT(disp.
Size() == Aop.
Height(),
"Sizes don't match!");
186 Aop.
Mult(disp, fullrhs);
189 for (
int k = 0; k < eliminators.Size(); ++k)
229 delete explicit_projector;
230 delete h_explicit_projector;
243 "Wrong number of dofs for elimination!");
245 for (
int i = 0; i < lagrange_dofs.
Size(); ++i)
247 lagrange_dofs[i] = i;
269 for (
int k = 0; k < constraint_rowstarts.
Size() - 1; ++k)
271 int constraint_size = constraint_rowstarts[k + 1] -
272 constraint_rowstarts[k];
278 for (
int i = constraint_rowstarts[k]; i < constraint_rowstarts[k + 1]; ++i)
280 lagrange_dofs[i - constraint_rowstarts[k]] = i;
281 for (
int jptr = I[i]; jptr < I[i + 1]; ++jptr)
285 if (std::abs(val) > 1.e-12 && secondary_dofs.
Find(j) == -1)
287 secondary_dofs[i - constraint_rowstarts[k]] = j;
293 for (
int i = constraint_rowstarts[k]; i < constraint_rowstarts[k + 1]; ++i)
295 MFEM_ASSERT(secondary_dofs[i - constraint_rowstarts[k]] >= 0,
296 "Secondary dofs don't match rows!");
297 for (
int jptr = I[i]; jptr < I[i + 1]; ++jptr)
300 if (secondary_dofs.
Find(j) == -1)
347 hA.
Mult(-1.0, rtilde, 1.0, temprhs);
384 MPI_Comm_rank(
A.GetComm(), &rank);
385 MPI_Comm_size(
A.GetComm(), &size);
387 int constraint_running_total = 0;
388 int local_constraints =
B.
Height();
389 MPI_Scan(&local_constraints, &constraint_running_total, 1, MPI_INT,
390 MPI_SUM,
A.GetComm());
391 int global_constraints = 0;
392 if (rank == size - 1) { global_constraints = constraint_running_total; }
393 MPI_Bcast(&global_constraints, 1, MPI_INT, size - 1,
A.GetComm());
397 HYPRE_BigInt row_starts[2] = { constraint_running_total - local_constraints,
398 constraint_running_total
402 row_starts, col_starts, &
B);
479 penalized_rhs += temp;
504class IdentitySolver :
public Solver
507 IdentitySolver(
int size) :
Solver(size) { }
508 void Mult(
const Vector& x, Vector& y)
const { y = x; }
509 void SetOperator(
const Operator& op) { }
512void SchurConstrainedSolver::Initialize()
521 tr_B =
new TransposeOperator(&
B);
535 primal_pc(&primal_pc_),
552 primal_pc(&primal_pc_),
601 if (
GetComm() != MPI_COMM_NULL)
641 h_primal_pc->SetPrintLevel(0);
644 h_primal_pc->SetSystemsOptions(
dimension, reorder);
658 schur_mat =
ParMult(scaledB, scaledBT);
662 h_dual_pc->SetPrintLevel(0);
678void ConstrainedSolver::Initialize()
718 for (
int i = 0; i <
f.Size(); ++i)
723 for (
int i = 0; i <
B.
Height(); ++i)
730 for (
int i = 0; i <
f.Size(); ++i)
734 for (
int i = 0; i <
B.
Height(); ++i)
741 Vector& x_and_lambda)
const
758 int node,
bool parallel,
int d=0)
767 const int vdof = pfespace->
DofToVDof(node, d);
772 MFEM_ABORT(
"Asked for parallel form of serial object!");
794 std::map<int, int> dof_bconstraint;
797 std::vector<std::map<int, int> > constraints;
798 int n_bconstraints = 0;
800 for (
int att : constrained_att)
803 std::set<int> constrained_tdofs;
804 for (
int i = 0; i < fespace.
GetNBE(); ++i)
817 if (tdof >= 0) { constrained_tdofs.insert(tdof); }
823 for (
auto k : constrained_tdofs)
825 auto it = dof_bconstraint.find(k);
826 if (it == dof_bconstraint.end())
829 dof_bconstraint[k] = n_bconstraints++;
830 constraints.emplace_back();
831 constraints.back()[att] = n_rows++;
836 constraints[it->second][att] = n_rows++;
844 std::map<int, int> reorder_rows;
847 constraint_rowstarts.
Append(0);
848 for (
auto& it : dof_bconstraint)
850 int bconstraint_index = it.second;
851 bool nconstraint =
false;
852 for (
auto& att_it : constraints[bconstraint_index])
854 auto rrit = reorder_rows.find(att_it.second);
855 if (rrit == reorder_rows.end())
858 reorder_rows[att_it.second] = new_row++;
861 if (nconstraint) { constraint_rowstarts.
Append(new_row); }
863 MFEM_VERIFY(new_row == n_rows,
"Remapping failed!");
864 for (
auto& constraint_map : constraints)
866 for (
auto& it : constraint_map)
868 it.second = reorder_rows[it.second];
878 std::map<int, int> node_visits;
879 for (
int i = 0; i < fespace.
GetNBE(); ++i)
890 MFEM_VERIFY(dofs.
Size() == nodes.
Size(),
891 "Something wrong in finite element space!");
893 for (
int j = 0; j < dofs.
Size(); ++j)
904 auto nv_it = node_visits.find(truek);
905 if (nv_it == node_visits.end())
907 node_visits[truek] = 1;
911 node_visits[truek]++;
913 int visits = node_visits[truek];
914 int bconstraint = dof_bconstraint[truek];
915 int row = constraints[bconstraint][att];
916 for (
int d = 0; d <
dim; ++d)
922 mout->
Add(row, inner_truek, nor[d]);
932 mout->
Set(row, inner_truek,
933 scaling * pv + (1.0 / visits) * nor[d]);
952 constraint_rowstarts,
true);
int FindSorted(const T &el) const
Do bisection search for 'el' in a sorted array; return -1 if not found.
void Sort()
Sorts the array in ascending order. This requires operator< to be defined for T.
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
int Size() const
Return the logical size of the array.
void DeleteAll()
Delete the whole array.
int Find(const T &el) const
Return the first index where 'el' is found; return -1 if not found.
int Append(const T &el)
Append element 'el' to array, resize if necessary.
void Unique()
Removes duplicities from a sorted array. This requires operator== to be defined for T.
T * HostReadWrite()
Shortcut for mfem::ReadWrite(a.GetMemory(), a.Size(), false).
A class to handle Block diagonal preconditioners in a matrix-free implementation.
void SetDiagonalBlock(int iblock, Operator *op)
Add a square block op in the block-entry (iblock, iblock).
Array< int > & RowOffsets()
Return the row offsets for block starts.
void SetBlock(int iRow, int iCol, Operator *op, real_t c=1.0)
Add a block op in the block-entry (iblock, jblock).
Operator & GetBlock(int i, int j)
Return a reference to block i,j.
An abstract class to solve the constrained system subject to the constraint .
virtual void LagrangeSystemMult(const Vector &f_and_r, Vector &x_and_lambda) const
Solve for (x, lambda) given (f, r)
void GetMultiplierSolution(Vector &lambda) const
Return the Lagrange multiplier solution in lambda.
virtual void SetConstraintRHS(const Vector &r)
Set the right-hand side r for the constraint B x = r.
ConstrainedSolver(MPI_Comm comm, Operator &A_, Operator &B_)
virtual void Mult(const Vector &f, Vector &x) const override
Solve for given .
Data type dense matrix using column-major storage.
void Mult(const real_t *x, real_t *y) const
Matrix vector multiplication.
void MultTranspose(const real_t *x, real_t *y) const
Multiply a vector with the transpose matrix.
void Transpose()
(*this) = (*this)^t
real_t * GetData() const
Returns the matrix data array.
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
real_t * HostReadWrite()
Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), false).
SparseMatrix * AssembleExact() const
Assemble this projector as a (processor-local) SparseMatrix.
void MultTranspose(const Vector &x, Vector &y) const override
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
void RecoverMultiplier(const Vector &primalrhs, const Vector &primalvars, Vector &lm) const
void Mult(const Vector &x, Vector &y) const override
Operator application: y=A(x).
EliminationProjection(const Operator &A, Array< Eliminator * > &eliminators)
void BuildGTilde(const Vector &g, Vector >ilde) const
void Mult(const Vector &x, Vector &y) const override
Solve for given .
Array< Eliminator * > eliminators
void BuildExplicitOperator()
Internal utility routine; assembles eliminated matrix explicitly.
EliminationSolver(HypreParMatrix &A, SparseMatrix &B, Array< int > &primary_dofs, Array< int > &secondary_dofs)
Constructor, with explicit splitting into primary/secondary dofs.
EliminationProjection * projector
virtual Solver * BuildPreconditioner() const =0
Build preconditioner for eliminated system.
virtual IterativeSolver * BuildKrylov() const =0
Select krylov solver for eliminated system.
HypreParMatrix * h_explicit_operator
Perform elimination of a single constraint.
void Eliminate(const Vector &in, Vector &out) const
void ExplicitAssembly(DenseMatrix &mat) const
Return explicitly assembled in mat.
void LagrangeSecondaryTranspose(const Vector &in, Vector &out) const
Transpose of LagrangeSecondary()
void LagrangeSecondary(const Vector &in, Vector &out) const
Maps Lagrange multipliers to secondary dofs, applies .
Eliminator(const SparseMatrix &B, const Array< int > &lagrange_dofs, const Array< int > &primary_tdofs, const Array< int > &secondary_tdofs)
const Array< int > & PrimaryDofs() const
const Array< int > & SecondaryDofs() const
void EliminateTranspose(const Vector &in, Vector &out) const
Transpose of Eliminate(), applies .
const Array< int > & LagrangeDofs() const
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th boundary fac...
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
int GetNBE() const
Returns number of boundary elements in the mesh.
int GetBdrAttribute(int i) const
ElementTransformation * GetBdrElementTransformation(int i) const
Returns ElementTransformation for the i-th boundary element.
DofTransformation * GetBdrElementDofs(int bel, Array< int > &dofs) const
Returns indices of degrees of freedom for boundary element 'bel'. The returned indices are offsets in...
int GetVDim() const
Returns vector dimension.
int DofToVDof(int dof, int vd, int ndofs=-1) const
Compute a single vdof corresponding to the index dof and the vector index vd.
Abstract class for all finite elements.
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
virtual void Mult(const Vector &b, Vector &x) const
Iterative solution of the linear system using the GMRES method.
The BoomerAMG solver in hypre.
Wrapper for hypre's ParCSR matrix class.
void GetDiag(Vector &diag) const
Get the local diagonal of the matrix.
void InvScaleRows(const Vector &s)
Scale the local row i by 1./s(i)
HYPRE_BigInt * GetRowStarts() const
Return the parallel row partitioning array.
HYPRE_BigInt GetGlobalNumRows() const
Return the global number of rows.
HYPRE_Int Mult(HypreParVector &x, HypreParVector &y, real_t alpha=1.0, real_t beta=0.0) const
Computes y = alpha * A * x + beta * y.
MPI_Comm GetComm() const
MPI communicator.
void EliminateZeroRows()
If a row contains only zeros, set its diagonal to 1.
HYPRE_BigInt * RowPart()
Returns the row partitioning.
HYPRE_BigInt M() const
Returns the global number of rows.
HypreParMatrix * Transpose() const
Returns the transpose of *this.
Class for an integration rule - an Array of IntegrationPoint.
Abstract base class for iterative solver.
real_t abs_tol
Absolute tolerance.
real_t rel_tol
Relative tolerance.
PrintLevel print_options
Output behavior for the iterative solver.
real_t GetFinalNorm() const
Returns the final residual norm after termination of the solver during the last call to Mult().
virtual void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
void SetRelTol(real_t rtol)
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
int GetNumIterations() const
Returns the number of iterations taken during the last call to Mult()
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
int max_iter
Limit for the number of iterations the solver is allowed to do.
void SetMaxIter(int max_it)
MPI_Comm GetComm() const
Return the associated MPI communicator, or MPI_COMM_NULL if no communicator is set.
real_t GetInitialNorm() const
Returns the initial residual norm from the last call to Mult().
bool GetConverged() const
Returns true if the last call to Mult() converged successfully.
void SetAbsTol(real_t atol)
virtual bool Factor(int m, real_t TOL=0.0)
Compute the LU factorization of the current matrix.
virtual void Solve(int m, int n, real_t *X) const
int width
Dimension of the input / number of columns in the matrix.
Operator(int s=0)
Construct a square Operator with given size s (default 0).
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
int height
Dimension of the output / number of rows in the matrix.
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Abstract parallel finite element space.
int GetLocalTDofNumber(int ldof) const
HypreParMatrix * penalized_mat
virtual IterativeSolver * BuildKrylov() const =0
Select krylov solver for penalized system.
void Mult(const Vector &x, Vector &y) const override
Solve for given .
virtual Solver * BuildPreconditioner() const =0
Build preconditioner for penalized system.
~PenaltyConstrainedSolver()
PenaltyConstrainedSolver(HypreParMatrix &A, SparseMatrix &B, real_t penalty_)
void Initialize(HypreParMatrix &A, HypreParMatrix &B, HypreParMatrix &D)
SchurConstrainedHypreSolver(MPI_Comm comm, HypreParMatrix &hA_, HypreParMatrix &hB_, Solver *prec=nullptr, int dimension=0, bool reorder=false)
virtual ~SchurConstrainedHypreSolver()
Solve constrained system by solving original mixed system; see ConstrainedSolver.
BlockDiagonalPreconditioner * block_pc
SchurConstrainedSolver(MPI_Comm comm, Operator &A_, Operator &B_, Solver &primal_pc_)
virtual ~SchurConstrainedSolver()
virtual void LagrangeSystemMult(const Vector &x, Vector &y) const override
Solve for (x, lambda) given (f, r)
virtual void SetOperator(const Operator &op)=0
Set/update the solver for the given operator.
void GetSubMatrix(const Array< int > &rows, const Array< int > &cols, DenseMatrix &subm) const
void Add(const int i, const int j, const real_t val)
void SetColPtr(const int row) const
Initialize the SparseMatrix for fast access to the entries of the given row which becomes the "curren...
real_t & SearchRow(const int col)
Perform a fast search for an entry in the "current row". See SetColPtr().
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
virtual void Finalize(int skip_zeros=1)
Finalize the matrix initialization, switching the storage format from LIL to CSR.
void Set(const int i, const int j, const real_t val)
void SetSubVector(const Array< int > &dofs, const real_t value)
Set the entries listed in dofs to the given value.
void AddElementVector(const Array< int > &dofs, const Vector &elemvect)
Add elements of the elemvect Vector to the entries listed in dofs. Negative dof values cause the -dof...
int Size() const
Returns the size of the vector.
void SetSize(int s)
Resize the vector to size s.
real_t * GetData() const
Return a pointer to the beginning of the Vector data.
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
constexpr int dimension
This example only works in 3D. Kernels for 2D are not implemented.
void CalcOrtho(const DenseMatrix &J, Vector &n)
int CanonicalNodeNumber(FiniteElementSpace &fespace, int node, bool parallel, int d=0)
SparseMatrix * BuildNormalConstraints(FiniteElementSpace &fespace, Array< int > &constrained_att, Array< int > &constraint_rowstarts, bool parallel)
Build a matrix constraining normal components to zero.
void RAP(const DenseMatrix &A, const DenseMatrix &P, DenseMatrix &RAP)
HypreParMatrix * ParAdd(const HypreParMatrix *A, const HypreParMatrix *B)
Returns the matrix A + B.
HypreParMatrix * ParMult(const HypreParMatrix *A, const HypreParMatrix *B, bool own_matrix)
SparseMatrix * ParBuildNormalConstraints(ParFiniteElementSpace &fespace, Array< int > &constrained_att, Array< int > &constraint_rowstarts)
Parallel wrapper for BuildNormalConstraints.
std::function< real_t(const Vector &)> f(real_t mass_coeff)