MFEM  v4.5.2
Finite element discretization library
constraints.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #include "constraints.hpp"
13 
14 #include "../fem/fespace.hpp"
15 #include "../fem/pfespace.hpp"
16 
17 #include <set>
18 
19 namespace mfem
20 {
21 
22 Eliminator::Eliminator(const SparseMatrix& B, const Array<int>& lagrange_tdofs_,
23  const Array<int>& primary_tdofs_,
24  const Array<int>& secondary_tdofs_)
25  :
26  lagrange_tdofs(lagrange_tdofs_),
27  primary_tdofs(primary_tdofs_),
28  secondary_tdofs(secondary_tdofs_)
29 {
30  MFEM_VERIFY(lagrange_tdofs.Size() == secondary_tdofs.Size(),
31  "Dof sizes don't match!");
32 
33  Bp.SetSize(lagrange_tdofs.Size(), primary_tdofs.Size());
34  B.GetSubMatrix(lagrange_tdofs, primary_tdofs, Bp);
35 
36  Bs.SetSize(lagrange_tdofs.Size(), secondary_tdofs.Size());
37  B.GetSubMatrix(lagrange_tdofs, secondary_tdofs, Bs);
38  BsT.Transpose(Bs);
39 
40  ipiv.SetSize(Bs.Height());
41  Bsinverse.data = Bs.HostReadWrite();
42  Bsinverse.ipiv = ipiv.HostReadWrite();
43  Bsinverse.Factor(Bs.Height());
44 
45  ipivT.SetSize(Bs.Height());
46  BsTinverse.data = BsT.HostReadWrite();
47  BsTinverse.ipiv = ipivT.HostReadWrite();
48  BsTinverse.Factor(Bs.Height());
49 }
50 
51 void Eliminator::Eliminate(const Vector& vin, Vector& vout) const
52 {
53  Bp.Mult(vin, vout);
54  Bsinverse.Solve(Bs.Height(), 1, vout.GetData());
55  vout *= -1.0;
56 }
57 
58 void Eliminator::EliminateTranspose(const Vector& vin, Vector& vout) const
59 {
60  Vector work(vin);
61  BsTinverse.Solve(Bs.Height(), 1, work.GetData());
62  Bp.MultTranspose(work, vout);
63  vout *= -1.0;
64 }
65 
66 void Eliminator::LagrangeSecondary(const Vector& vin, Vector& vout) const
67 {
68  vout = vin;
69  Bsinverse.Solve(Bs.Height(), 1, vout.GetData());
70 }
71 
73  Vector& vout) const
74 {
75  vout = vin;
76  BsTinverse.Solve(Bs.Height(), 1, vout.GetData());
77 }
78 
80 {
81  mat.SetSize(Bp.Height(), Bp.Width());
82  mat = Bp;
83  Bsinverse.Solve(Bs.Height(), Bp.Width(), mat.GetData());
84  mat *= -1.0;
85 }
86 
88  Array<Eliminator*>& eliminators_)
89  :
90  Operator(A.Height()),
91  Aop(A),
92  eliminators(eliminators_)
93 {
94 }
95 
96 void EliminationProjection::Mult(const Vector& vin, Vector& vout) const
97 {
98  MFEM_ASSERT(vin.Size() == width, "Wrong vector size!");
99  MFEM_ASSERT(vout.Size() == height, "Wrong vector size!");
100 
101  vout = vin;
102 
103  for (int k = 0; k < eliminators.Size(); ++k)
104  {
105  Eliminator* elim = eliminators[k];
106  Vector subvec_in;
107  Vector subvec_out(elim->SecondaryDofs().Size());
108  vin.GetSubVector(elim->PrimaryDofs(), subvec_in);
109  elim->Eliminate(subvec_in, subvec_out);
110  vout.SetSubVector(elim->SecondaryDofs(), subvec_out);
111  }
112 }
113 
115 {
116  MFEM_ASSERT(vin.Size() == height, "Wrong vector size!");
117  MFEM_ASSERT(vout.Size() == width, "Wrong vector size!");
118 
119  vout = vin;
120 
121  for (int k = 0; k < eliminators.Size(); ++k)
122  {
123  Eliminator* elim = eliminators[k];
124  Vector subvec_in;
125  Vector subvec_out(elim->PrimaryDofs().Size());
126  vin.GetSubVector(elim->SecondaryDofs(), subvec_in);
127  elim->EliminateTranspose(subvec_in, subvec_out);
128  vout.AddElementVector(elim->PrimaryDofs(), subvec_out);
129  vout.SetSubVector(elim->SecondaryDofs(), 0.0);
130  }
131 }
132 
134 {
135  SparseMatrix * mat = new SparseMatrix(height, width);
136 
137  for (int i = 0; i < height; ++i)
138  {
139  mat->Add(i, i, 1.0);
140  }
141 
142  for (int k = 0; k < eliminators.Size(); ++k)
143  {
144  Eliminator* elim = eliminators[k];
145  DenseMatrix mat_k;
146  elim->ExplicitAssembly(mat_k);
147  for (int iz = 0; iz < elim->SecondaryDofs().Size(); ++iz)
148  {
149  int i = elim->SecondaryDofs()[iz];
150  for (int jz = 0; jz < elim->PrimaryDofs().Size(); ++jz)
151  {
152  int j = elim->PrimaryDofs()[jz];
153  mat->Add(i, j, mat_k(iz, jz));
154  }
155  mat->Set(i, i, 0.0);
156  }
157  }
158 
159  mat->Finalize();
160  return mat;
161 }
162 
163 void EliminationProjection::BuildGTilde(const Vector& r, Vector& rtilde) const
164 {
165  MFEM_ASSERT(rtilde.Size() == Aop.Height(), "Sizes don't match!");
166 
167  rtilde = 0.0;
168  for (int k = 0; k < eliminators.Size(); ++k)
169  {
170  Eliminator* elim = eliminators[k];
171  Vector subr;
172  r.GetSubVector(elim->LagrangeDofs(), subr);
173  Vector bsinvr(subr.Size());
174  elim->LagrangeSecondary(subr, bsinvr);
175  rtilde.AddElementVector(elim->SecondaryDofs(), bsinvr);
176  }
177 }
178 
180  const Vector& disprhs, const Vector& disp, Vector& lagrangem) const
181 {
182  lagrangem = 0.0;
183  MFEM_ASSERT(disp.Size() == Aop.Height(), "Sizes don't match!");
184 
185  Vector fullrhs(Aop.Height());
186  Aop.Mult(disp, fullrhs);
187  fullrhs -= disprhs;
188  fullrhs *= -1.0;
189  for (int k = 0; k < eliminators.Size(); ++k)
190  {
191  Eliminator* elim = eliminators[k];
192  Vector localsec;
193  fullrhs.GetSubVector(elim->SecondaryDofs(), localsec);
194  Vector locallagrange(localsec.Size());
195  elim->LagrangeSecondaryTranspose(localsec, locallagrange);
196  lagrangem.AddElementVector(elim->LagrangeDofs(), locallagrange);
197  }
198 }
199 
200 #ifdef MFEM_USE_MPI
201 
203 {
204  delete h_explicit_operator;
205  for (auto elim : eliminators)
206  {
207  delete elim;
208  }
209  delete projector;
210  delete prec;
211  delete krylov;
212 }
213 
215 {
216  SparseMatrix * explicit_projector = projector->AssembleExact();
217  HypreParMatrix * h_explicit_projector =
219  hA.GetRowStarts(), explicit_projector);
220  h_explicit_projector->CopyRowStarts();
221  h_explicit_projector->CopyColStarts();
222 
223  h_explicit_operator = RAP(&hA, h_explicit_projector);
224  // next line because of square projector
228 
229  delete explicit_projector;
230  delete h_explicit_projector;
231 }
232 
234  Array<int>& primary_dofs,
235  Array<int>& secondary_dofs)
236  :
237  ConstrainedSolver(A.GetComm(), A, B),
238  hA(A),
239  krylov(nullptr),
240  prec(nullptr)
241 {
242  MFEM_VERIFY(secondary_dofs.Size() == B.Height(),
243  "Wrong number of dofs for elimination!");
244  Array<int> lagrange_dofs(secondary_dofs.Size());
245  for (int i = 0; i < lagrange_dofs.Size(); ++i)
246  {
247  lagrange_dofs[i] = i;
248  }
249  eliminators.Append(new Eliminator(B, lagrange_dofs, primary_dofs,
250  secondary_dofs));
253 }
254 
256  Array<int>& constraint_rowstarts)
257  :
258  ConstrainedSolver(A.GetComm(), A, B),
259  hA(A),
260  krylov(nullptr),
261  prec(nullptr)
262 {
263  if (!B.Empty())
264  {
265  int * I = B.GetI();
266  int * J = B.GetJ();
267  double * data = B.GetData();
268 
269  for (int k = 0; k < constraint_rowstarts.Size() - 1; ++k)
270  {
271  int constraint_size = constraint_rowstarts[k + 1] -
272  constraint_rowstarts[k];
273  Array<int> lagrange_dofs(constraint_size);
274  Array<int> primary_dofs;
275  Array<int> secondary_dofs(constraint_size);
276  secondary_dofs = -1;
277  // loop through rows, identify one secondary dof for each row
278  for (int i = constraint_rowstarts[k]; i < constraint_rowstarts[k + 1]; ++i)
279  {
280  lagrange_dofs[i - constraint_rowstarts[k]] = i;
281  for (int jptr = I[i]; jptr < I[i + 1]; ++jptr)
282  {
283  int j = J[jptr];
284  double val = data[jptr];
285  if (std::abs(val) > 1.e-12 && secondary_dofs.Find(j) == -1)
286  {
287  secondary_dofs[i - constraint_rowstarts[k]] = j;
288  break;
289  }
290  }
291  }
292  // loop through rows again, assigning non-secondary dofs as primary
293  for (int i = constraint_rowstarts[k]; i < constraint_rowstarts[k + 1]; ++i)
294  {
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)
298  {
299  int j = J[jptr];
300  if (secondary_dofs.Find(j) == -1)
301  {
302  primary_dofs.Append(j);
303  }
304  }
305  }
306  primary_dofs.Sort();
307  primary_dofs.Unique();
308  eliminators.Append(new Eliminator(B, lagrange_dofs, primary_dofs,
309  secondary_dofs));
310  }
311  }
314 }
315 
316 void EliminationSolver::Mult(const Vector& rhs, Vector& sol) const
317 {
318  if (!prec)
319  {
321  }
322  else
323  {
325  }
326  if (!krylov)
327  {
328  krylov = BuildKrylov();
331  }
336 
337  Vector rtilde(rhs.Size());
338  if (constraint_rhs.Size() > 0)
339  {
341  }
342  else
343  {
344  rtilde = 0.0;
345  }
346  Vector temprhs(rhs);
347  hA.Mult(-1.0, rtilde, 1.0, temprhs);
348 
349  Vector reducedrhs(rhs.Size());
350  projector->MultTranspose(temprhs, reducedrhs);
351  Vector reducedsol(rhs.Size());
352  reducedsol = 0.0;
353  krylov->Mult(reducedrhs, reducedsol);
358 
359  projector->Mult(reducedsol, sol);
361  sol += rtilde;
362 }
363 
364 void PenaltyConstrainedSolver::Initialize(HypreParMatrix& A, HypreParMatrix& B,
365  HypreParMatrix& D)
366 {
367  HypreParMatrix * hBTB = RAP(&D, &B);
368  // this matrix doesn't get cleanly deleted?
369  // (hypre comm pkg)
370  penalized_mat = ParAdd(&A, hBTB);
371  delete hBTB;
372 }
373 
375  HypreParMatrix& A, SparseMatrix& B, double penalty_)
376  :
377  ConstrainedSolver(A.GetComm(), A, B),
378  penalty(B.Height()),
379  constraintB(B),
380  krylov(nullptr),
381  prec(nullptr)
382 {
383  int rank, size;
384  MPI_Comm_rank(A.GetComm(), &rank);
385  MPI_Comm_size(A.GetComm(), &size);
386 
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());
394 
395  HYPRE_BigInt glob_num_rows = global_constraints;
396  HYPRE_BigInt glob_num_cols = A.N();
397  HYPRE_BigInt row_starts[2] = { constraint_running_total - local_constraints,
398  constraint_running_total
399  };
400  HYPRE_BigInt col_starts[2] = { A.ColPart()[0], A.ColPart()[1] };
401  HypreParMatrix hB(A.GetComm(), glob_num_rows, glob_num_cols,
402  row_starts, col_starts, &B);
403  hB.CopyRowStarts();
404  hB.CopyColStarts();
405  penalty=penalty_;
407  HypreParMatrix hD(hB.GetComm(), hB.M(), hB.RowPart(), &D);
408  hD.CopyRowStarts();
409  hD.CopyColStarts();
410  Initialize(A, hB, hD);
411 }
412 
414  HypreParMatrix& A, HypreParMatrix& B, double penalty_)
415  :
416  ConstrainedSolver(A.GetComm(), A, B),
417  penalty(B.Height()),
418  constraintB(B),
419  krylov(nullptr),
420  prec(nullptr)
421 {
422  penalty=penalty_;
424  HypreParMatrix hD(B.GetComm(), B.M(), B.RowPart(), &D);
425  hD.CopyRowStarts();
426  hD.CopyColStarts();
427  Initialize(A, B, hD);
428 }
429 
431  HypreParMatrix& A, HypreParMatrix& B, Vector& penalty_)
432  :
433  ConstrainedSolver(A.GetComm(), A, B),
434  penalty(penalty_),
435  constraintB(B),
436  krylov(nullptr),
437  prec(nullptr)
438 {
439  SparseMatrix D(penalty_);
440  HypreParMatrix hD(B.GetComm(), B.M(), B.RowPart(), &D);
441  hD.CopyRowStarts();
442  hD.CopyColStarts();
443  Initialize(A, B, hD);
444 }
445 
447 {
448  delete penalized_mat;
449  delete prec;
450  delete krylov;
451 }
452 
454 {
455  if (!prec)
456  {
458  }
459  else
460  {
462  }
463  if (!krylov)
464  {
465  krylov = BuildKrylov();
468  }
469 
470  // form penalized right-hand side
471  Vector penalized_rhs(b);
472  if (constraint_rhs.Size() > 0)
473  {
474  Vector temp_rhs(constraint_rhs.Size());
476  D.Mult(constraint_rhs, temp_rhs);
477  Vector temp(x.Size());
478  constraintB.MultTranspose(temp_rhs, temp);
479  penalized_rhs += temp;
480  }
481 
482  // actually solve
487  krylov->Mult(penalized_rhs, x);
492 
494  if (constraint_rhs.Size() > 0)
495  {
497  }
499 }
500 
501 #endif
502 
503 /// because IdentityOperator isn't a Solver
504 class IdentitySolver : public Solver
505 {
506 public:
507  IdentitySolver(int size) : Solver(size) { }
508  void Mult(const Vector& x, Vector& y) const { y = x; }
509  void SetOperator(const Operator& op) { }
510 };
511 
512 void SchurConstrainedSolver::Initialize()
513 {
514  offsets[0] = 0;
515  offsets[1] = A.Height();
516  offsets[2] = A.Height() + B.Height();
517 
518  block_op = new BlockOperator(offsets);
519  block_op->SetBlock(0, 0, &A);
520  block_op->SetBlock(1, 0, &B);
521  tr_B = new TransposeOperator(&B);
522  block_op->SetBlock(0, 1, tr_B);
523 
524  block_pc = new BlockDiagonalPreconditioner(block_op->RowOffsets()),
525  rel_tol = 1.e-6;
526 }
527 
528 #ifdef MFEM_USE_MPI
530  Operator& A_, Operator& B_,
531  Solver& primal_pc_)
532  :
533  ConstrainedSolver(comm, A_, B_),
534  offsets(3),
535  primal_pc(&primal_pc_),
536  dual_pc(nullptr)
537 {
538  Initialize();
540  dual_pc = new IdentitySolver(block_op->RowOffsets()[2] -
541  block_op->RowOffsets()[1]);
544 }
545 #endif
546 
548  Solver& primal_pc_)
549  :
550  ConstrainedSolver(A_, B_),
551  offsets(3),
552  primal_pc(&primal_pc_),
553  dual_pc(nullptr)
554 {
555  Initialize();
557  dual_pc = new IdentitySolver(block_op->RowOffsets()[2] -
558  block_op->RowOffsets()[1]);
561 }
562 
563 #ifdef MFEM_USE_MPI
564 // protected constructor
566  Operator& B_)
567  :
568  ConstrainedSolver(comm, A_, B_),
569  offsets(3),
570  primal_pc(nullptr),
571  dual_pc(nullptr)
572 {
573  Initialize();
574 }
575 #endif
576 
577 // protected constructor
579  :
580  ConstrainedSolver(A_, B_),
581  offsets(3),
582  primal_pc(nullptr),
583  dual_pc(nullptr)
584 {
585  Initialize();
586 }
587 
589 {
590  delete block_op;
591  delete tr_B;
592  delete block_pc;
593  delete dual_pc;
594 }
595 
597  Vector& y) const
598 {
599  GMRESSolver * gmres;
600 #ifdef MFEM_USE_MPI
601  if (GetComm() != MPI_COMM_NULL)
602  {
603  gmres = new GMRESSolver(GetComm());
604  }
605  else
606 #endif
607  {
608  gmres = new GMRESSolver;
609  }
610  gmres->SetOperator(*block_op);
611  gmres->SetRelTol(rel_tol);
612  gmres->SetAbsTol(abs_tol);
613  gmres->SetMaxIter(max_iter);
615  gmres->SetPreconditioner(
616  const_cast<BlockDiagonalPreconditioner&>(*block_pc));
617 
618  gmres->Mult(x, y);
619  final_iter = gmres->GetNumIterations();
620  converged = gmres->GetConverged();
621  initial_norm = gmres->GetInitialNorm();
622  final_norm = gmres->GetFinalNorm();
623  delete gmres;
624 }
625 
626 #ifdef MFEM_USE_MPI
628  HypreParMatrix& hA_,
629  HypreParMatrix& hB_,
630  Solver * prec,
631  int dimension,
632  bool reorder)
633  :
634  SchurConstrainedSolver(comm, hA_, hB_),
635  hA(hA_),
636  hB(hB_)
637 {
638  if (prec == nullptr)
639  {
640  auto h_primal_pc = new HypreBoomerAMG(hA);
641  h_primal_pc->SetPrintLevel(0);
642  if (dimension > 0)
643  {
644  h_primal_pc->SetSystemsOptions(dimension, reorder);
645  }
646  primal_pc = h_primal_pc;
647  }
648  else
649  {
650  primal_pc = prec;
651  }
652 
653  HypreParMatrix * scaledB = new HypreParMatrix(hB);
654  Vector diagA;
655  hA.GetDiag(diagA);
656  HypreParMatrix * scaledBT = scaledB->Transpose();
657  scaledBT->InvScaleRows(diagA);
658  schur_mat = ParMult(scaledB, scaledBT);
659  schur_mat->CopyRowStarts();
660  schur_mat->CopyColStarts();
661  auto h_dual_pc = new HypreBoomerAMG(*schur_mat);
662  h_dual_pc->SetPrintLevel(0);
663  dual_pc = h_dual_pc;
664  delete scaledB;
665  delete scaledBT;
666 
669 }
670 
672 {
673  delete schur_mat;
674  delete primal_pc;
675 }
676 #endif
677 
678 void ConstrainedSolver::Initialize()
679 {
680  height = A.Height() + B.Height();
681  width = A.Width() + B.Height();
682 
683  workb.SetSize(A.Height());
684  workx.SetSize(A.Height());
686  constraint_rhs = 0.0;
688 }
689 
690 #ifdef MFEM_USE_MPI
692  :
693  IterativeSolver(comm), A(A_), B(B_)
694 {
695  Initialize();
696 }
697 #endif
698 
700  :
701  A(A_), B(B_)
702 {
703  Initialize();
704 }
705 
707 {
708  MFEM_VERIFY(r.Size() == multiplier_sol.Size(), "Vector is wrong size!");
709  constraint_rhs = r;
710 }
711 
712 void ConstrainedSolver::Mult(const Vector& f, Vector &x) const
713 {
714  Vector pworkb(A.Height() + B.Height());
715  Vector pworkx(A.Height() + B.Height());
716  pworkb = 0.0;
717  pworkx = 0.0;
718  for (int i = 0; i < f.Size(); ++i)
719  {
720  pworkb(i) = f(i);
721  pworkx(i) = x(i);
722  }
723  for (int i = 0; i < B.Height(); ++i)
724  {
725  pworkb(f.Size() + i) = constraint_rhs(i);
726  }
727 
728  LagrangeSystemMult(pworkb, pworkx);
729 
730  for (int i = 0; i < f.Size(); ++i)
731  {
732  x(i) = pworkx(i);
733  }
734  for (int i = 0; i < B.Height(); ++i)
735  {
736  multiplier_sol(i) = pworkx(f.Size() + i);
737  }
738 }
739 
741  Vector& x_and_lambda) const
742 {
743  workb.MakeRef(const_cast<Vector&>(f_and_r), 0);
744  workx.MakeRef(x_and_lambda, 0);
745  Vector ref_constraint_rhs(f_and_r.GetData() + A.Height(), B.Height());
746  constraint_rhs = ref_constraint_rhs;
747  Mult(workb, workx);
748  Vector ref_constraint_sol(x_and_lambda.GetData() + A.Height(), B.Height());
749  GetMultiplierSolution(ref_constraint_sol);
750 }
751 
752 /* Helper routine to reduce code duplication - given a node (which MFEM
753  sometimes calls a "dof"), this returns what normal people call a dof but
754  which MFEM sometimes calls a "vdof" - note that MFEM's naming conventions
755  regarding this are not entirely consistent. In parallel, this always
756  returns the "truedof" in parallel numbering. */
758  int node, bool parallel, int d=0)
759 {
760 #ifdef MFEM_USE_MPI
761  if (parallel)
762  {
763  ParFiniteElementSpace* pfespace =
764  dynamic_cast<ParFiniteElementSpace*>(&fespace);
765  if (pfespace)
766  {
767  const int vdof = pfespace->DofToVDof(node, d);
768  return pfespace->GetLocalTDofNumber(vdof);
769  }
770  else
771  {
772  MFEM_ABORT("Asked for parallel form of serial object!");
773  return -1;
774  }
775  }
776  else
777 #endif
778  {
779  return fespace.DofToVDof(node, d);
780  }
781 }
782 
784  Array<int>& constrained_att,
785  Array<int>& constraint_rowstarts,
786  bool parallel)
787 {
788  int dim = fespace.GetVDim();
789 
790  // dof_constraint maps a dof (column of the constraint matrix) to
791  // a block-constraint
792  // the indexing is by tdof, but a single tdof uniquely identifies a node
793  // so we only store one tdof independent of dimension
794  std::map<int, int> dof_bconstraint;
795  // constraints[j] is a map from attribute to row number,
796  // the j itself is the index of a block-constraint
797  std::vector<std::map<int, int> > constraints;
798  int n_bconstraints = 0;
799  int n_rows = 0;
800  for (int att : constrained_att)
801  {
802  // identify tdofs on constrained boundary
803  std::set<int> constrained_tdofs;
804  for (int i = 0; i < fespace.GetNBE(); ++i)
805  {
806  if (fespace.GetBdrAttribute(i) == att)
807  {
808  Array<int> nodes;
809  // get nodes on boundary (MFEM sometimes calls these dofs, what
810  // we call dofs it calls vdofs)
811  fespace.GetBdrElementDofs(i, nodes);
812  for (auto k : nodes)
813  {
814  // get the (local) dof number corresponding to
815  // the x-coordinate dof for node k
816  int tdof = CanonicalNodeNumber(fespace, k, parallel);
817  if (tdof >= 0) { constrained_tdofs.insert(tdof); }
818  }
819  }
820  }
821  // fill in the maps identifying which constraints (rows) correspond to
822  // which tdofs
823  for (auto k : constrained_tdofs)
824  {
825  auto it = dof_bconstraint.find(k);
826  if (it == dof_bconstraint.end())
827  {
828  // build new block constraint
829  dof_bconstraint[k] = n_bconstraints++;
830  constraints.emplace_back();
831  constraints.back()[att] = n_rows++;
832  }
833  else
834  {
835  // add tdof to existing block constraint
836  constraints[it->second][att] = n_rows++;
837  }
838  }
839  }
840 
841  // reorder so block-constraints eliminated together are grouped together in
842  // adjacent rows
843  {
844  std::map<int, int> reorder_rows;
845  int new_row = 0;
846  constraint_rowstarts.DeleteAll();
847  constraint_rowstarts.Append(0);
848  for (auto& it : dof_bconstraint)
849  {
850  int bconstraint_index = it.second;
851  bool nconstraint = false;
852  for (auto& att_it : constraints[bconstraint_index])
853  {
854  auto rrit = reorder_rows.find(att_it.second);
855  if (rrit == reorder_rows.end())
856  {
857  nconstraint = true;
858  reorder_rows[att_it.second] = new_row++;
859  }
860  }
861  if (nconstraint) { constraint_rowstarts.Append(new_row); }
862  }
863  MFEM_VERIFY(new_row == n_rows, "Remapping failed!");
864  for (auto& constraint_map : constraints)
865  {
866  for (auto& it : constraint_map)
867  {
868  it.second = reorder_rows[it.second];
869  }
870  }
871  }
872 
873  SparseMatrix * mout = new SparseMatrix(n_rows, fespace.GetTrueVSize());
874 
875  // fill in constraint matrix with normal vector information
876  Vector nor(dim);
877  // how many times we have seen a node (key is truek)
878  std::map<int, int> node_visits;
879  for (int i = 0; i < fespace.GetNBE(); ++i)
880  {
881  int att = fespace.GetBdrAttribute(i);
882  if (constrained_att.FindSorted(att) != -1)
883  {
885  const FiniteElement * fe = fespace.GetBE(i);
886  const IntegrationRule& nodes = fe->GetNodes();
887 
888  Array<int> dofs;
889  fespace.GetBdrElementDofs(i, dofs);
890  MFEM_VERIFY(dofs.Size() == nodes.Size(),
891  "Something wrong in finite element space!");
892 
893  for (int j = 0; j < dofs.Size(); ++j)
894  {
895  Tr->SetIntPoint(&nodes[j]);
896  // the normal returned in the next line is scaled by h, which is
897  // probably what we want in most applications
898  CalcOrtho(Tr->Jacobian(), nor);
899 
900  int k = dofs[j];
901  int truek = CanonicalNodeNumber(fespace, k, parallel);
902  if (truek >= 0)
903  {
904  auto nv_it = node_visits.find(truek);
905  if (nv_it == node_visits.end())
906  {
907  node_visits[truek] = 1;
908  }
909  else
910  {
911  node_visits[truek]++;
912  }
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)
917  {
918  int inner_truek = CanonicalNodeNumber(fespace, k,
919  parallel, d);
920  if (visits == 1)
921  {
922  mout->Add(row, inner_truek, nor[d]);
923  }
924  else
925  {
926  mout->SetColPtr(row);
927  const double pv = mout->SearchRow(inner_truek);
928  const double scaling = ((double) (visits - 1)) /
929  ((double) visits);
930  // incremental average, based on how many times
931  // this node has been visited
932  mout->Set(row, inner_truek,
933  scaling * pv + (1.0 / visits) * nor[d]);
934  }
935 
936  }
937  }
938  }
939  }
940  }
941  mout->Finalize();
942 
943  return mout;
944 }
945 
946 #ifdef MFEM_USE_MPI
948  Array<int>& constrained_att,
949  Array<int>& constraint_rowstarts)
950 {
951  return BuildNormalConstraints(fespace, constrained_att,
952  constraint_rowstarts, true);
953 }
954 #endif
955 
956 }
Abstract class for all finite elements.
Definition: fe_base.hpp:232
HypreParMatrix & hA
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 ...
Definition: operator.hpp:93
T * HostReadWrite()
Shortcut for mfem::ReadWrite(a.GetMemory(), a.Size(), false).
Definition: array.hpp:327
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
void Mult(const Vector &x, Vector &y) const override
Solve for given .
void Unique()
Removes duplicities from a sorted array. This requires operator== to be defined for T...
Definition: array.hpp:259
Array< int > & RowOffsets()
Return the row offsets for block starts.
SparseMatrix * ParBuildNormalConstraints(ParFiniteElementSpace &fespace, Array< int > &constrained_att, Array< int > &constraint_rowstarts)
Parallel wrapper for BuildNormalConstraints.
Solve constrained system by solving original mixed system; see ConstrainedSolver. ...
EliminationSolver(HypreParMatrix &A, SparseMatrix &B, Array< int > &primary_dofs, Array< int > &secondary_dofs)
Constructor, with explicit splitting into primary/secondary dofs.
void Mult(const Vector &x, Vector &y) const override
Operator application: y=A(x).
Definition: constraints.cpp:96
void Initialize(HypreParMatrix &A, HypreParMatrix &B, HypreParMatrix &D)
double GetInitialNorm() const
Returns the initial residual norm from the last call to Mult().
Definition: solvers.hpp:258
void RecoverMultiplier(const Vector &primalrhs, const Vector &primalvars, Vector &lm) const
virtual void Finalize(int skip_zeros=1)
Finalize the matrix initialization, switching the storage format from LIL to CSR. ...
Definition: sparsemat.hpp:546
SparseMatrix * BuildNormalConstraints(FiniteElementSpace &fespace, Array< int > &constrained_att, Array< int > &constraint_rowstarts, bool parallel)
Build a matrix constraining normal components to zero.
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:512
double & SearchRow(const int col)
Perform a fast search for an entry in the "current row". See SetColPtr().
Definition: sparsemat.hpp:846
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
SchurConstrainedHypreSolver(MPI_Comm comm, HypreParMatrix &hA_, HypreParMatrix &hB_, Solver *prec=nullptr, int dimension=0, bool reorder=false)
int CanonicalNodeNumber(FiniteElementSpace &fespace, int node, bool parallel, int d=0)
void SetColPtr(const int row) const
Initialize the SparseMatrix for fast access to the entries of the given row which becomes the "curren...
Definition: sparsemat.hpp:792
An abstract class to solve the constrained system subject to the constraint .
Definition: constraints.hpp:53
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
virtual DofTransformation * GetBdrElementDofs(int bel, Array< int > &dofs) const
Returns indices of degrees of freedom for boundary element &#39;bel&#39;.
Definition: fespace.cpp:2814
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:93
virtual void Solve(int m, int n, double *X) const
Definition: densemat.cpp:3420
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:974
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
virtual void LagrangeSystemMult(const Vector &f_and_r, Vector &x_and_lambda) const
Solve for (x, lambda) given (f, r)
void ExplicitAssembly(DenseMatrix &mat) const
Return explicitly assembled in mat.
Definition: constraints.cpp:79
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
Abstract parallel finite element space.
Definition: pfespace.hpp:28
double * HostReadWrite()
Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), false).
Definition: densemat.hpp:486
const Array< int > & SecondaryDofs() const
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition: vector.cpp:548
HypreParMatrix * h_explicit_operator
Eliminator(const SparseMatrix &B, const Array< int > &lagrange_dofs, const Array< int > &primary_tdofs, const Array< int > &secondary_tdofs)
Definition: constraints.cpp:22
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 GetSubMatrix(const Array< int > &rows, const Array< int > &cols, DenseMatrix &subm) const
Definition: sparsemat.cpp:2886
void CalcOrtho(const DenseMatrix &J, Vector &n)
Definition: densemat.cpp:2693
HypreParMatrix * ParMult(const HypreParMatrix *A, const HypreParMatrix *B, bool own_matrix)
Definition: hypre.cpp:2850
virtual Solver * BuildPreconditioner() const =0
Build preconditioner for penalized system.
void DeleteAll()
Delete the whole array.
Definition: array.hpp:851
Array< Eliminator * > eliminators
Operator & GetBlock(int i, int j)
Return a reference to block i,j.
double GetFinalNorm() const
Returns the final residual norm after termination of the solver during the last call to Mult()...
Definition: solvers.hpp:265
ConstrainedSolver(MPI_Comm comm, Operator &A_, Operator &B_)
SparseMatrix * AssembleExact() const
Assemble this projector as a (processor-local) SparseMatrix.
virtual void Mult(const Vector &f, Vector &x) const override
Solve for given .
void EliminateZeroRows()
If a row contains only zeros, set its diagonal to 1.
Definition: hypre.hpp:799
virtual void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:179
void Set(const int i, const int j, const double val)
Definition: sparsemat.cpp:2768
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1590
void EliminateTranspose(const Vector &in, Vector &out) const
Transpose of Eliminate(), applies .
Definition: constraints.cpp:58
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition: solvers.cpp:71
double f(const Vector &xvec)
Definition: lor_mms.hpp:32
double * GetData() const
Returns the matrix data array.
Definition: densemat.hpp:115
Data type sparse matrix.
Definition: sparsemat.hpp:50
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:756
Perform elimination of a single constraint.
virtual IterativeSolver * BuildKrylov() const =0
Select krylov solver for penalized system.
ElementTransformation * GetBdrElementTransformation(int i) const
Returns ElementTransformation for the i-th boundary element.
Definition: fespace.hpp:656
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
Definition: densemat.cpp:215
int GetNumIterations() const
Returns the number of iterations taken during the last call to Mult()
Definition: solvers.hpp:250
double b
Definition: lissajous.cpp:42
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint, using the method SetIntPoint().
Definition: eltrans.hpp:119
int Find(const T &el) const
Return the first index where &#39;el&#39; is found; return -1 if not found.
Definition: array.hpp:815
void SetMaxIter(int max_it)
Definition: solvers.hpp:201
virtual void SetConstraintRHS(const Vector &r)
Set the right-hand side r for the constraint B x = r.
void Sort()
Sorts the array in ascending order. This requires operator< to be defined for T.
Definition: array.hpp:251
double * data
Definition: densemat.hpp:622
void GetDiag(Vector &diag) const
Get the local diagonal of the matrix.
Definition: hypre.cpp:1481
void RAP(const DenseMatrix &A, const DenseMatrix &P, DenseMatrix &RAP)
Definition: densemat.cpp:3252
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:173
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...
Definition: vector.cpp:640
void BuildExplicitOperator()
Internal utility routine; assembles eliminated matrix explicitly.
void Add(const int i, const int j, const double val)
Definition: sparsemat.cpp:2787
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:208
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:581
MPI_Comm GetComm() const
MPI communicator.
Definition: hypre.hpp:534
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:590
void LagrangeSecondary(const Vector &in, Vector &out) const
Maps Lagrange multipliers to secondary dofs, applies .
Definition: constraints.cpp:66
void SetAbsTol(double atol)
Definition: solvers.hpp:200
PrintLevel print_options
Output behavior for the iterative solver.
Definition: solvers.hpp:137
void SetRelTol(double rtol)
Definition: solvers.hpp:199
PenaltyConstrainedSolver(HypreParMatrix &A, SparseMatrix &B, double penalty_)
Abstract base class for iterative solver.
Definition: solvers.hpp:66
void Transpose()
(*this) = (*this)^t
Definition: densemat.cpp:1420
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
EliminationProjection(const Operator &A, Array< Eliminator *> &eliminators)
Definition: constraints.cpp:87
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:684
bool GetConverged() const
Returns true if the last call to Mult() converged successfully.
Definition: solvers.hpp:252
int max_iter
Limit for the number of iterations the solver is allowed to do.
Definition: solvers.hpp:151
HYPRE_Int HYPRE_BigInt
int GetLocalTDofNumber(int ldof) const
Definition: pfespace.cpp:1062
GMRES method.
Definition: solvers.hpp:525
SchurConstrainedSolver(MPI_Comm comm, Operator &A_, Operator &B_, Solver &primal_pc_)
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
HypreParMatrix * ParAdd(const HypreParMatrix *A, const HypreParMatrix *B)
Returns the matrix A + B.
Definition: hypre.cpp:2810
virtual Solver * BuildPreconditioner() const =0
Build preconditioner for eliminated system.
HypreParMatrix * Transpose() const
Returns the transpose of *this.
Definition: hypre.cpp:1608
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
EliminationProjection * projector
const Array< int > & LagrangeDofs() const
double rel_tol
Relative tolerance.
Definition: solvers.hpp:154
void Mult(const Vector &x, Vector &y) const override
Solve for given .
void InvScaleRows(const Vector &s)
Scale the local row i by 1./s(i)
Definition: hypre.cpp:2053
int GetBdrAttribute(int i) const
Definition: fespace.hpp:661
IterativeSolver * krylov
int dim
Definition: ex24.cpp:53
void BuildGTilde(const Vector &g, Vector &gtilde) const
virtual bool Factor(int m, double TOL=0.0)
Compute the LU factorization of the current matrix.
Definition: densemat.cpp:3269
int GetNBE() const
Returns number of boundary elements in the mesh.
Definition: fespace.hpp:623
constexpr int dimension
This example only works in 3D. Kernels for 2D are not implemented.
Definition: hooke.cpp:45
virtual void LagrangeSystemMult(const Vector &x, Vector &y) const override
Solve for (x, lambda) given (f, r)
void Eliminate(const Vector &in, Vector &out) const
Definition: constraints.cpp:51
virtual void SetOperator(const Operator &op)=0
Set/update the solver for the given operator.
virtual IterativeSolver * BuildKrylov() const =0
Select krylov solver for eliminated system.
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
MPI_Comm GetComm() const
Return the associated MPI communicator, or MPI_COMM_NULL if no communicator is set.
Definition: solvers.hpp:292
int DofToVDof(int dof, int vd, int ndofs=-1) const
Definition: fespace.cpp:251
Vector data type.
Definition: vector.hpp:60
const Array< int > & PrimaryDofs() const
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
Definition: vector.hpp:576
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:173
BlockDiagonalPreconditioner * block_pc
Base class for solvers.
Definition: operator.hpp:682
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition: fe_base.hpp:388
void GetMultiplierSolution(Vector &lambda) const
Return the Lagrange multiplier solution in lambda.
Definition: constraints.hpp:76
TransposeOperator * tr_B
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th boundary fac...
Definition: fespace.cpp:3102
HYPRE_BigInt GetGlobalNumRows() const
Return the global number of rows.
Definition: hypre.hpp:635
HYPRE_Int Mult(HypreParVector &x, HypreParVector &y, double alpha=1.0, double beta=0.0) const
Computes y = alpha * A * x + beta * y.
Definition: hypre.cpp:1733
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28
void SetBlock(int iRow, int iCol, Operator *op, double c=1.0)
Add a block op in the block-entry (iblock, jblock).
double abs_tol
Absolute tolerance.
Definition: solvers.hpp:157
void LagrangeSecondaryTranspose(const Vector &in, Vector &out) const
Transpose of LagrangeSecondary()
Definition: constraints.cpp:72
void SetDiagonalBlock(int iblock, Operator *op)
Add a square block op in the block-entry (iblock, iblock).
HYPRE_BigInt * GetRowStarts() const
Return the parallel row partitioning array.
Definition: hypre.hpp:645