MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
constraints.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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
19namespace mfem
20{
21
22Eliminator::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
51void 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
58void 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
66void 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
96void 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{
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
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 real_t * 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 real_t 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
316void EliminationSolver::Mult(const Vector& rhs, Vector& sol) const
317{
318 if (!prec)
319 {
321 }
322 else
323 {
325 }
326 if (!krylov)
327 {
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
364void PenaltyConstrainedSolver::Initialize(HypreParMatrix& A, HypreParMatrix& B,
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, real_t 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, real_t 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 {
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
504class IdentitySolver : public Solver
505{
506public:
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
512void 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(
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
678void ConstrainedSolver::Initialize()
679{
680 height = A.Height() + B.Height();
681 width = A.Width() + B.Height();
682
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
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 real_t pv = mout->SearchRow(inner_truek);
928 const real_t scaling = ((real_t) (visits - 1)) /
929 ((real_t) 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}
int FindSorted(const T &el) const
Do bisection search for 'el' in a sorted array; return -1 if not found.
Definition array.hpp:838
void Sort()
Sorts the array in ascending order. This requires operator< to be defined for T.
Definition array.hpp:261
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition array.hpp:697
int Size() const
Return the logical size of the array.
Definition array.hpp:144
void DeleteAll()
Delete the whole array.
Definition array.hpp:864
int Find(const T &el) const
Return the first index where 'el' is found; return -1 if not found.
Definition array.hpp:828
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Definition array.hpp:769
void Unique()
Removes duplicities from a sorted array. This requires operator== to be defined for T.
Definition array.hpp:269
T * HostReadWrite()
Shortcut for mfem::ReadWrite(a.GetMemory(), a.Size(), false).
Definition array.hpp:337
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.
Definition densemat.hpp:24
void Mult(const real_t *x, real_t *y) const
Matrix vector multiplication.
Definition densemat.cpp:220
void MultTranspose(const real_t *x, real_t *y) const
Multiply a vector with the transpose matrix.
Definition densemat.cpp:262
void Transpose()
(*this) = (*this)^t
real_t * GetData() const
Returns the matrix data array.
Definition densemat.hpp:115
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition densemat.hpp:105
real_t * HostReadWrite()
Shortcut for mfem::ReadWrite(GetMemory(), TotalSize(), false).
Definition densemat.hpp:486
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint,...
Definition eltrans.hpp:119
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition eltrans.hpp:93
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 &gtilde) const
void Mult(const Vector &x, Vector &y) const override
Solve for given .
IterativeSolver * krylov
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
real_t * data
Definition densemat.hpp:629
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th boundary fac...
Definition fespace.cpp:3204
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition fespace.hpp:716
int GetNBE() const
Returns number of boundary elements in the mesh.
Definition fespace.hpp:749
int GetBdrAttribute(int i) const
Definition fespace.hpp:787
ElementTransformation * GetBdrElementTransformation(int i) const
Returns ElementTransformation for the i-th boundary element.
Definition fespace.hpp:782
DofTransformation * GetBdrElementDofs(int bel, Array< int > &dofs) const
Returns indices of degrees of freedom for boundary element 'bel'. The returned indices are offsets in...
Definition fespace.cpp:2950
int GetVDim() const
Returns vector dimension.
Definition fespace.hpp:706
int DofToVDof(int dof, int vd, int ndofs=-1) const
Compute a single vdof corresponding to the index dof and the vector index vd.
Definition fespace.cpp:249
Abstract class for all finite elements.
Definition fe_base.hpp:239
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition fe_base.hpp:395
GMRES method.
Definition solvers.hpp:547
virtual void Mult(const Vector &b, Vector &x) const
Iterative solution of the linear system using the GMRES method.
Definition solvers.cpp:980
The BoomerAMG solver in hypre.
Definition hypre.hpp:1691
Wrapper for hypre's ParCSR matrix class.
Definition hypre.hpp:388
void GetDiag(Vector &diag) const
Get the local diagonal of the matrix.
Definition hypre.cpp:1557
void InvScaleRows(const Vector &s)
Scale the local row i by 1./s(i)
Definition hypre.cpp:2135
HYPRE_BigInt * GetRowStarts() const
Return the parallel row partitioning array.
Definition hypre.hpp:689
HYPRE_BigInt GetGlobalNumRows() const
Return the global number of rows.
Definition hypre.hpp:679
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.
Definition hypre.cpp:1815
MPI_Comm GetComm() const
MPI communicator.
Definition hypre.hpp:578
void EliminateZeroRows()
If a row contains only zeros, set its diagonal to 1.
Definition hypre.hpp:843
HYPRE_BigInt * RowPart()
Returns the row partitioning.
Definition hypre.hpp:613
HYPRE_BigInt M() const
Returns the global number of rows.
Definition hypre.hpp:627
HypreParMatrix * Transpose() const
Returns the transpose of *this.
Definition hypre.cpp:1684
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
Abstract base class for iterative solver.
Definition solvers.hpp:67
real_t abs_tol
Absolute tolerance.
Definition solvers.hpp:158
real_t rel_tol
Relative tolerance.
Definition solvers.hpp:155
PrintLevel print_options
Output behavior for the iterative solver.
Definition solvers.hpp:138
real_t GetFinalNorm() const
Returns the final residual norm after termination of the solver during the last call to Mult().
Definition solvers.hpp:275
virtual void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
Definition solvers.cpp:179
void SetRelTol(real_t rtol)
Definition solvers.hpp:209
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition solvers.cpp:173
int GetNumIterations() const
Returns the number of iterations taken during the last call to Mult()
Definition solvers.hpp:260
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition solvers.cpp:71
int max_iter
Limit for the number of iterations the solver is allowed to do.
Definition solvers.hpp:152
void SetMaxIter(int max_it)
Definition solvers.hpp:211
MPI_Comm GetComm() const
Return the associated MPI communicator, or MPI_COMM_NULL if no communicator is set.
Definition solvers.hpp:302
real_t GetInitialNorm() const
Returns the initial residual norm from the last call to Mult().
Definition solvers.hpp:268
bool GetConverged() const
Returns true if the last call to Mult() converged successfully.
Definition solvers.hpp:262
void SetAbsTol(real_t atol)
Definition solvers.hpp:210
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
Abstract operator.
Definition operator.hpp:25
int width
Dimension of the input / number of columns in the matrix.
Definition operator.hpp:28
Operator(int s=0)
Construct a square Operator with given size s (default 0).
Definition operator.hpp:59
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition operator.hpp:66
int height
Dimension of the output / number of rows in the matrix.
Definition operator.hpp:27
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().
Definition operator.hpp:72
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
Abstract parallel finite element space.
Definition pfespace.hpp:29
int GetLocalTDofNumber(int ldof) const
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(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)
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 void LagrangeSystemMult(const Vector &x, Vector &y) const override
Solve for (x, lambda) given (f, r)
Base class for solvers.
Definition operator.hpp:683
virtual void SetOperator(const Operator &op)=0
Set/update the solver for the given operator.
Data type sparse matrix.
Definition sparsemat.hpp:51
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)
Vector data type.
Definition vector.hpp:80
void SetSubVector(const Array< int > &dofs, const real_t value)
Set the entries listed in dofs to the given value.
Definition vector.cpp:604
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:670
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
real_t * GetData() const
Return a pointer to the beginning of the Vector data.
Definition vector.hpp:227
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition vector.cpp:578
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
Definition vector.hpp:602
int dim
Definition ex24.cpp:53
constexpr int dimension
This example only works in 3D. Kernels for 2D are not implemented.
Definition hooke.cpp:45
HYPRE_Int HYPRE_BigInt
real_t b
Definition lissajous.cpp:42
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.
Definition hypre.cpp:2909
HypreParMatrix * ParMult(const HypreParMatrix *A, const HypreParMatrix *B, bool own_matrix)
Definition hypre.cpp:2949
SparseMatrix * ParBuildNormalConstraints(ParFiniteElementSpace &fespace, Array< int > &constrained_att, Array< int > &constraint_rowstarts)
Parallel wrapper for BuildNormalConstraints.
float real_t
Definition config.hpp:43
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30