1 // Copyright (c) 2010-2021, 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
Eliminator(const SparseMatrix & B,const Array<int> & lagrange_tdofs_,const Array<int> & primary_tdofs_,const Array<int> & secondary_tdofs_)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
Eliminate(const Vector & in,Vector & out) const51 void Eliminator::Eliminate(const Vector& in, Vector& out) const
52 {
53 Bp.Mult(in, out);
54 Bsinverse.Solve(Bs.Height(), 1, out);
55 out *= -1.0;
56 }
57
EliminateTranspose(const Vector & in,Vector & out) const58 void Eliminator::EliminateTranspose(const Vector& in, Vector& out) const
59 {
60 Vector work(in);
61 BsTinverse.Solve(Bs.Height(), 1, work);
62 Bp.MultTranspose(work, out);
63 out *= -1.0;
64 }
65
LagrangeSecondary(const Vector & in,Vector & out) const66 void Eliminator::LagrangeSecondary(const Vector& in, Vector& out) const
67 {
68 out = in;
69 Bsinverse.Solve(Bs.Height(), 1, out);
70 }
71
LagrangeSecondaryTranspose(const Vector & in,Vector & out) const72 void Eliminator::LagrangeSecondaryTranspose(const Vector& in, Vector& out) const
73 {
74 out = in;
75 BsTinverse.Solve(Bs.Height(), 1, out);
76 }
77
ExplicitAssembly(DenseMatrix & mat) const78 void Eliminator::ExplicitAssembly(DenseMatrix& mat) const
79 {
80 mat.SetSize(Bp.Height(), Bp.Width());
81 mat = Bp;
82 Bsinverse.Solve(Bs.Height(), Bp.Width(), mat.GetData());
83 mat *= -1.0;
84 }
85
EliminationProjection(const Operator & A,Array<Eliminator * > & eliminators_)86 EliminationProjection::EliminationProjection(const Operator& A,
87 Array<Eliminator*>& eliminators_)
88 :
89 Operator(A.Height()),
90 Aop(A),
91 eliminators(eliminators_)
92 {
93 }
94
Mult(const Vector & in,Vector & out) const95 void EliminationProjection::Mult(const Vector& in, Vector& out) const
96 {
97 MFEM_ASSERT(in.Size() == width, "Wrong vector size!");
98 MFEM_ASSERT(out.Size() == height, "Wrong vector size!");
99
100 out = in;
101
102 for (int k = 0; k < eliminators.Size(); ++k)
103 {
104 Eliminator* elim = eliminators[k];
105 Vector subvec_in;
106 Vector subvec_out(elim->SecondaryDofs().Size());
107 in.GetSubVector(elim->PrimaryDofs(), subvec_in);
108 elim->Eliminate(subvec_in, subvec_out);
109 out.SetSubVector(elim->SecondaryDofs(), subvec_out);
110 }
111 }
112
MultTranspose(const Vector & in,Vector & out) const113 void EliminationProjection::MultTranspose(const Vector& in, Vector& out) const
114 {
115 MFEM_ASSERT(in.Size() == height, "Wrong vector size!");
116 MFEM_ASSERT(out.Size() == width, "Wrong vector size!");
117
118 out = in;
119
120 for (int k = 0; k < eliminators.Size(); ++k)
121 {
122 Eliminator* elim = eliminators[k];
123 Vector subvec_in;
124 Vector subvec_out(elim->PrimaryDofs().Size());
125 in.GetSubVector(elim->SecondaryDofs(), subvec_in);
126 elim->EliminateTranspose(subvec_in, subvec_out);
127 out.AddElementVector(elim->PrimaryDofs(), subvec_out);
128 out.SetSubVector(elim->SecondaryDofs(), 0.0);
129 }
130 }
131
AssembleExact() const132 SparseMatrix * EliminationProjection::AssembleExact() const
133 {
134 SparseMatrix * out = new SparseMatrix(height, width);
135
136 for (int i = 0; i < height; ++i)
137 {
138 out->Add(i, i, 1.0);
139 }
140
141 for (int k = 0; k < eliminators.Size(); ++k)
142 {
143 Eliminator* elim = eliminators[k];
144 DenseMatrix mat;
145 elim->ExplicitAssembly(mat);
146 for (int iz = 0; iz < elim->SecondaryDofs().Size(); ++iz)
147 {
148 int i = elim->SecondaryDofs()[iz];
149 for (int jz = 0; jz < elim->PrimaryDofs().Size(); ++jz)
150 {
151 int j = elim->PrimaryDofs()[jz];
152 out->Add(i, j, mat(iz, jz));
153 }
154 out->Set(i, i, 0.0);
155 }
156 }
157
158 out->Finalize();
159 return out;
160 }
161
BuildGTilde(const Vector & r,Vector & rtilde) const162 void EliminationProjection::BuildGTilde(const Vector& r, Vector& rtilde) const
163 {
164 MFEM_ASSERT(rtilde.Size() == Aop.Height(), "Sizes don't match!");
165
166 rtilde = 0.0;
167 for (int k = 0; k < eliminators.Size(); ++k)
168 {
169 Eliminator* elim = eliminators[k];
170 Vector subr;
171 r.GetSubVector(elim->LagrangeDofs(), subr);
172 Vector bsinvr(subr.Size());
173 elim->LagrangeSecondary(subr, bsinvr);
174 rtilde.AddElementVector(elim->SecondaryDofs(), bsinvr);
175 }
176 }
177
RecoverMultiplier(const Vector & disprhs,const Vector & disp,Vector & lagrangem) const178 void EliminationProjection::RecoverMultiplier(
179 const Vector& disprhs, const Vector& disp, Vector& lagrangem) const
180 {
181 lagrangem = 0.0;
182 MFEM_ASSERT(disp.Size() == Aop.Height(), "Sizes don't match!");
183
184 Vector fullrhs(Aop.Height());
185 Aop.Mult(disp, fullrhs);
186 fullrhs -= disprhs;
187 fullrhs *= -1.0;
188 for (int k = 0; k < eliminators.Size(); ++k)
189 {
190 Eliminator* elim = eliminators[k];
191 Vector localsec;
192 fullrhs.GetSubVector(elim->SecondaryDofs(), localsec);
193 Vector locallagrange(localsec.Size());
194 elim->LagrangeSecondaryTranspose(localsec, locallagrange);
195 lagrangem.AddElementVector(elim->LagrangeDofs(), locallagrange);
196 }
197 }
198
199 #ifdef MFEM_USE_MPI
200
~EliminationSolver()201 EliminationSolver::~EliminationSolver()
202 {
203 delete h_explicit_operator;
204 for (auto elim : eliminators)
205 {
206 delete elim;
207 }
208 delete projector;
209 delete prec;
210 delete krylov;
211 }
212
BuildExplicitOperator()213 void EliminationSolver::BuildExplicitOperator()
214 {
215 SparseMatrix * explicit_projector = projector->AssembleExact();
216 HypreParMatrix * h_explicit_projector =
217 new HypreParMatrix(hA.GetComm(), hA.GetGlobalNumRows(),
218 hA.GetRowStarts(), explicit_projector);
219 h_explicit_projector->CopyRowStarts();
220 h_explicit_projector->CopyColStarts();
221
222 h_explicit_operator = RAP(&hA, h_explicit_projector);
223 // next line because of square projector
224 h_explicit_operator->EliminateZeroRows();
225 h_explicit_operator->CopyRowStarts();
226 h_explicit_operator->CopyColStarts();
227
228 delete explicit_projector;
229 delete h_explicit_projector;
230 }
231
EliminationSolver(HypreParMatrix & A,SparseMatrix & B,Array<int> & primary_dofs,Array<int> & secondary_dofs)232 EliminationSolver::EliminationSolver(HypreParMatrix& A, SparseMatrix& B,
233 Array<int>& primary_dofs,
234 Array<int>& secondary_dofs)
235 :
236 ConstrainedSolver(A.GetComm(), A, B),
237 hA(A),
238 krylov(nullptr),
239 prec(nullptr)
240 {
241 MFEM_VERIFY(secondary_dofs.Size() == B.Height(),
242 "Wrong number of dofs for elimination!");
243 Array<int> lagrange_dofs(secondary_dofs.Size());
244 for (int i = 0; i < lagrange_dofs.Size(); ++i)
245 {
246 lagrange_dofs[i] = i;
247 }
248 eliminators.Append(new Eliminator(B, lagrange_dofs, primary_dofs,
249 secondary_dofs));
250 projector = new EliminationProjection(hA, eliminators);
251 BuildExplicitOperator();
252 }
253
EliminationSolver(HypreParMatrix & A,SparseMatrix & B,Array<int> & constraint_rowstarts)254 EliminationSolver::EliminationSolver(HypreParMatrix& A, SparseMatrix& B,
255 Array<int>& constraint_rowstarts)
256 :
257 ConstrainedSolver(A.GetComm(), A, B),
258 hA(A),
259 krylov(nullptr),
260 prec(nullptr)
261 {
262 if (!B.Empty())
263 {
264 int * I = B.GetI();
265 int * J = B.GetJ();
266 double * data = B.GetData();
267
268 for (int k = 0; k < constraint_rowstarts.Size() - 1; ++k)
269 {
270 int constraint_size = constraint_rowstarts[k + 1] -
271 constraint_rowstarts[k];
272 Array<int> lagrange_dofs(constraint_size);
273 Array<int> primary_dofs;
274 Array<int> secondary_dofs(constraint_size);
275 secondary_dofs = -1;
276 // loop through rows, identify one secondary dof for each row
277 for (int i = constraint_rowstarts[k]; i < constraint_rowstarts[k + 1]; ++i)
278 {
279 lagrange_dofs[i - constraint_rowstarts[k]] = i;
280 for (int jptr = I[i]; jptr < I[i + 1]; ++jptr)
281 {
282 int j = J[jptr];
283 double val = data[jptr];
284 if (std::abs(val) > 1.e-12 && secondary_dofs.Find(j) == -1)
285 {
286 secondary_dofs[i - constraint_rowstarts[k]] = j;
287 break;
288 }
289 }
290 }
291 // loop through rows again, assigning non-secondary dofs as primary
292 for (int i = constraint_rowstarts[k]; i < constraint_rowstarts[k + 1]; ++i)
293 {
294 MFEM_ASSERT(secondary_dofs[i - constraint_rowstarts[k]] >= 0,
295 "Secondary dofs don't match rows!");
296 for (int jptr = I[i]; jptr < I[i + 1]; ++jptr)
297 {
298 int j = J[jptr];
299 if (secondary_dofs.Find(j) == -1)
300 {
301 primary_dofs.Append(j);
302 }
303 }
304 }
305 primary_dofs.Sort();
306 primary_dofs.Unique();
307 eliminators.Append(new Eliminator(B, lagrange_dofs, primary_dofs,
308 secondary_dofs));
309 }
310 }
311 projector = new EliminationProjection(hA, eliminators);
312 BuildExplicitOperator();
313 }
314
Mult(const Vector & rhs,Vector & sol) const315 void EliminationSolver::Mult(const Vector& rhs, Vector& sol) const
316 {
317 if (!prec)
318 {
319 prec = BuildPreconditioner();
320 }
321 if (!krylov)
322 {
323 krylov = BuildKrylov();
324 krylov->SetOperator(*h_explicit_operator);
325 krylov->SetPreconditioner(*prec);
326 }
327 krylov->SetMaxIter(max_iter);
328 krylov->SetRelTol(rel_tol);
329 krylov->SetAbsTol(abs_tol);
330 krylov->SetPrintLevel(print_level);
331
332 Vector rtilde(rhs.Size());
333 if (constraint_rhs.Size() > 0)
334 {
335 projector->BuildGTilde(constraint_rhs, rtilde);
336 }
337 else
338 {
339 rtilde = 0.0;
340 }
341 Vector temprhs(rhs);
342 hA.Mult(-1.0, rtilde, 1.0, temprhs);
343
344 Vector reducedrhs(rhs.Size());
345 projector->MultTranspose(temprhs, reducedrhs);
346 Vector reducedsol(rhs.Size());
347 reducedsol = 0.0;
348 krylov->Mult(reducedrhs, reducedsol);
349 final_iter = krylov->GetNumIterations();
350 final_norm = krylov->GetFinalNorm();
351 converged = krylov->GetConverged();
352
353 projector->Mult(reducedsol, sol);
354 projector->RecoverMultiplier(temprhs, sol, multiplier_sol);
355 sol += rtilde;
356 }
357
Initialize(HypreParMatrix & A,HypreParMatrix & B)358 void PenaltyConstrainedSolver::Initialize(HypreParMatrix& A, HypreParMatrix& B)
359 {
360 HypreParMatrix * hBT = B.Transpose();
361 HypreParMatrix * hBTB = ParMult(hBT, &B, true);
362 // this matrix doesn't get cleanly deleted?
363 // (hypre comm pkg)
364 (*hBTB) *= penalty;
365 penalized_mat = ParAdd(&A, hBTB);
366 delete hBTB;
367 delete hBT;
368 }
369
PenaltyConstrainedSolver(HypreParMatrix & A,SparseMatrix & B,double penalty_)370 PenaltyConstrainedSolver::PenaltyConstrainedSolver(
371 HypreParMatrix& A, SparseMatrix& B, double penalty_)
372 :
373 ConstrainedSolver(A.GetComm(), A, B),
374 penalty(penalty_),
375 constraintB(B),
376 krylov(nullptr),
377 prec(nullptr)
378 {
379 int rank, size;
380 MPI_Comm_rank(A.GetComm(), &rank);
381 MPI_Comm_size(A.GetComm(), &size);
382
383 int constraint_running_total = 0;
384 int local_constraints = B.Height();
385 MPI_Scan(&local_constraints, &constraint_running_total, 1, MPI_INT,
386 MPI_SUM, A.GetComm());
387 int global_constraints = 0;
388 if (rank == size - 1) { global_constraints = constraint_running_total; }
389 MPI_Bcast(&global_constraints, 1, MPI_INT, size - 1, A.GetComm());
390
391 HYPRE_BigInt glob_num_rows = global_constraints;
392 HYPRE_BigInt glob_num_cols = A.N();
393 HYPRE_BigInt row_starts[2] = { constraint_running_total - local_constraints,
394 constraint_running_total
395 };
396 HYPRE_BigInt col_starts[2] = { A.ColPart()[0], A.ColPart()[1] };
397 HypreParMatrix hB(A.GetComm(), glob_num_rows, glob_num_cols,
398 row_starts, col_starts, &B);
399 hB.CopyRowStarts();
400 hB.CopyColStarts();
401 Initialize(A, hB);
402 }
403
PenaltyConstrainedSolver(HypreParMatrix & A,HypreParMatrix & B,double penalty_)404 PenaltyConstrainedSolver::PenaltyConstrainedSolver(
405 HypreParMatrix& A, HypreParMatrix& B, double penalty_)
406 :
407 ConstrainedSolver(A.GetComm(), A, B),
408 penalty(penalty_),
409 constraintB(B),
410 krylov(nullptr),
411 prec(nullptr)
412 {
413 Initialize(A, B);
414 }
415
~PenaltyConstrainedSolver()416 PenaltyConstrainedSolver::~PenaltyConstrainedSolver()
417 {
418 delete penalized_mat;
419 delete prec;
420 delete krylov;
421 }
422
Mult(const Vector & b,Vector & x) const423 void PenaltyConstrainedSolver::Mult(const Vector& b, Vector& x) const
424 {
425 if (!prec)
426 {
427 prec = BuildPreconditioner();
428 }
429 if (!krylov)
430 {
431 krylov = BuildKrylov();
432 krylov->SetOperator(*penalized_mat);
433 krylov->SetPreconditioner(*prec);
434 }
435
436 // form penalized right-hand side
437 Vector penalized_rhs(b);
438 if (constraint_rhs.Size() > 0)
439 {
440 Vector temp(x.Size());
441 constraintB.MultTranspose(constraint_rhs, temp);
442 temp *= penalty;
443 penalized_rhs += temp;
444 }
445
446 // actually solve
447 krylov->SetRelTol(rel_tol);
448 krylov->SetAbsTol(abs_tol);
449 krylov->SetMaxIter(max_iter);
450 krylov->SetPrintLevel(print_level);
451 krylov->Mult(penalized_rhs, x);
452 final_iter = krylov->GetNumIterations();
453 final_norm = krylov->GetFinalNorm();
454 converged = krylov->GetConverged();
455
456 constraintB.Mult(x, multiplier_sol);
457 if (constraint_rhs.Size() > 0)
458 {
459 multiplier_sol -= constraint_rhs;
460 }
461 multiplier_sol *= penalty;
462 }
463
464 #endif
465
466 /// because IdentityOperator isn't a Solver
467 class IdentitySolver : public Solver
468 {
469 public:
IdentitySolver(int size)470 IdentitySolver(int size) : Solver(size) { }
Mult(const Vector & x,Vector & y) const471 void Mult(const Vector& x, Vector& y) const { y = x; }
SetOperator(const Operator & op)472 void SetOperator(const Operator& op) { }
473 };
474
Initialize()475 void SchurConstrainedSolver::Initialize()
476 {
477 offsets[0] = 0;
478 offsets[1] = A.Height();
479 offsets[2] = A.Height() + B.Height();
480
481 block_op = new BlockOperator(offsets);
482 block_op->SetBlock(0, 0, &A);
483 block_op->SetBlock(1, 0, &B);
484 tr_B = new TransposeOperator(&B);
485 block_op->SetBlock(0, 1, tr_B);
486
487 block_pc = new BlockDiagonalPreconditioner(block_op->RowOffsets()),
488 rel_tol = 1.e-6;
489 }
490
491 #ifdef MFEM_USE_MPI
SchurConstrainedSolver(MPI_Comm comm,Operator & A_,Operator & B_,Solver & primal_pc_)492 SchurConstrainedSolver::SchurConstrainedSolver(MPI_Comm comm,
493 Operator& A_, Operator& B_,
494 Solver& primal_pc_)
495 :
496 ConstrainedSolver(comm, A_, B_),
497 offsets(3),
498 primal_pc(&primal_pc_),
499 dual_pc(nullptr)
500 {
501 Initialize();
502 primal_pc->SetOperator(block_op->GetBlock(0, 0));
503 dual_pc = new IdentitySolver(block_op->RowOffsets()[2] -
504 block_op->RowOffsets()[1]);
505 block_pc->SetDiagonalBlock(0, primal_pc);
506 block_pc->SetDiagonalBlock(1, dual_pc);
507 }
508 #endif
509
SchurConstrainedSolver(Operator & A_,Operator & B_,Solver & primal_pc_)510 SchurConstrainedSolver::SchurConstrainedSolver(Operator& A_, Operator& B_,
511 Solver& primal_pc_)
512 :
513 ConstrainedSolver(A_, B_),
514 offsets(3),
515 primal_pc(&primal_pc_),
516 dual_pc(nullptr)
517 {
518 Initialize();
519 primal_pc->SetOperator(block_op->GetBlock(0, 0));
520 dual_pc = new IdentitySolver(block_op->RowOffsets()[2] -
521 block_op->RowOffsets()[1]);
522 block_pc->SetDiagonalBlock(0, primal_pc);
523 block_pc->SetDiagonalBlock(1, dual_pc);
524 }
525
526 #ifdef MFEM_USE_MPI
527 // protected constructor
SchurConstrainedSolver(MPI_Comm comm,Operator & A_,Operator & B_)528 SchurConstrainedSolver::SchurConstrainedSolver(MPI_Comm comm, Operator& A_,
529 Operator& B_)
530 :
531 ConstrainedSolver(comm, A_, B_),
532 offsets(3),
533 primal_pc(nullptr),
534 dual_pc(nullptr)
535 {
536 Initialize();
537 }
538 #endif
539
540 // protected constructor
SchurConstrainedSolver(Operator & A_,Operator & B_)541 SchurConstrainedSolver::SchurConstrainedSolver(Operator& A_, Operator& B_)
542 :
543 ConstrainedSolver(A_, B_),
544 offsets(3),
545 primal_pc(nullptr),
546 dual_pc(nullptr)
547 {
548 Initialize();
549 }
550
~SchurConstrainedSolver()551 SchurConstrainedSolver::~SchurConstrainedSolver()
552 {
553 delete block_op;
554 delete tr_B;
555 delete block_pc;
556 delete dual_pc;
557 }
558
LagrangeSystemMult(const Vector & x,Vector & y) const559 void SchurConstrainedSolver::LagrangeSystemMult(const Vector& x,
560 Vector& y) const
561 {
562 GMRESSolver * gmres;
563 #ifdef MFEM_USE_MPI
564 if (GetComm() != MPI_COMM_NULL)
565 {
566 gmres = new GMRESSolver(GetComm());
567 }
568 else
569 #endif
570 {
571 gmres = new GMRESSolver;
572 }
573 gmres->SetOperator(*block_op);
574 gmres->SetRelTol(rel_tol);
575 gmres->SetAbsTol(abs_tol);
576 gmres->SetMaxIter(max_iter);
577 gmres->SetPrintLevel(print_level);
578 gmres->SetPreconditioner(
579 const_cast<BlockDiagonalPreconditioner&>(*block_pc));
580
581 gmres->Mult(x, y);
582 final_iter = gmres->GetNumIterations();
583 delete gmres;
584 }
585
586 #ifdef MFEM_USE_MPI
SchurConstrainedHypreSolver(MPI_Comm comm,HypreParMatrix & hA_,HypreParMatrix & hB_,int dimension,bool reorder)587 SchurConstrainedHypreSolver::SchurConstrainedHypreSolver(MPI_Comm comm,
588 HypreParMatrix& hA_,
589 HypreParMatrix& hB_,
590 int dimension,
591 bool reorder)
592 :
593 SchurConstrainedSolver(comm, hA_, hB_),
594 hA(hA_),
595 hB(hB_)
596 {
597 auto h_primal_pc = new HypreBoomerAMG(hA);
598 h_primal_pc->SetPrintLevel(0);
599 if (dimension > 0)
600 {
601 h_primal_pc->SetSystemsOptions(dimension, reorder);
602 }
603 primal_pc = h_primal_pc;
604
605 HypreParMatrix * scaledB = new HypreParMatrix(hB);
606 Vector diagA;
607 hA.GetDiag(diagA);
608 HypreParMatrix * scaledBT = scaledB->Transpose();
609 scaledBT->InvScaleRows(diagA);
610 schur_mat = ParMult(scaledB, scaledBT);
611 schur_mat->CopyRowStarts();
612 schur_mat->CopyColStarts();
613 auto h_dual_pc = new HypreBoomerAMG(*schur_mat);
614 h_dual_pc->SetPrintLevel(0);
615 dual_pc = h_dual_pc;
616 delete scaledB;
617 delete scaledBT;
618
619 block_pc->SetDiagonalBlock(0, primal_pc);
620 block_pc->SetDiagonalBlock(1, dual_pc);
621 }
622
~SchurConstrainedHypreSolver()623 SchurConstrainedHypreSolver::~SchurConstrainedHypreSolver()
624 {
625 delete schur_mat;
626 delete primal_pc;
627 }
628 #endif
629
Initialize()630 void ConstrainedSolver::Initialize()
631 {
632 height = A.Height() + B.Height();
633 width = A.Width() + B.Height();
634
635 workb.SetSize(A.Height());
636 workx.SetSize(A.Height());
637 constraint_rhs.SetSize(B.Height());
638 constraint_rhs = 0.0;
639 multiplier_sol.SetSize(B.Height());
640 }
641
642 #ifdef MFEM_USE_MPI
ConstrainedSolver(MPI_Comm comm,Operator & A_,Operator & B_)643 ConstrainedSolver::ConstrainedSolver(MPI_Comm comm, Operator& A_, Operator& B_)
644 :
645 IterativeSolver(comm), A(A_), B(B_)
646 {
647 Initialize();
648 }
649 #endif
650
ConstrainedSolver(Operator & A_,Operator & B_)651 ConstrainedSolver::ConstrainedSolver(Operator& A_, Operator& B_)
652 :
653 A(A_), B(B_)
654 {
655 Initialize();
656 }
657
SetConstraintRHS(const Vector & r)658 void ConstrainedSolver::SetConstraintRHS(const Vector& r)
659 {
660 MFEM_VERIFY(r.Size() == multiplier_sol.Size(), "Vector is wrong size!");
661 constraint_rhs = r;
662 }
663
Mult(const Vector & f,Vector & x) const664 void ConstrainedSolver::Mult(const Vector& f, Vector &x) const
665 {
666 Vector pworkb(A.Height() + B.Height());
667 Vector pworkx(A.Height() + B.Height());
668 pworkb = 0.0;
669 pworkx = 0.0;
670 for (int i = 0; i < f.Size(); ++i)
671 {
672 pworkb(i) = f(i);
673 pworkx(i) = x(i);
674 }
675 for (int i = 0; i < B.Height(); ++i)
676 {
677 pworkb(f.Size() + i) = constraint_rhs(i);
678 }
679
680 LagrangeSystemMult(pworkb, pworkx);
681
682 for (int i = 0; i < f.Size(); ++i)
683 {
684 x(i) = pworkx(i);
685 }
686 for (int i = 0; i < B.Height(); ++i)
687 {
688 multiplier_sol(i) = pworkx(f.Size() + i);
689 }
690 }
691
LagrangeSystemMult(const Vector & f_and_r,Vector & x_and_lambda) const692 void ConstrainedSolver::LagrangeSystemMult(const Vector& f_and_r,
693 Vector& x_and_lambda) const
694 {
695 workb.MakeRef(const_cast<Vector&>(f_and_r), 0);
696 workx.MakeRef(x_and_lambda, 0);
697 Vector ref_constraint_rhs(f_and_r.GetData() + A.Height(), B.Height());
698 constraint_rhs = ref_constraint_rhs;
699 Mult(workb, workx);
700 Vector ref_constraint_sol(x_and_lambda.GetData() + A.Height(), B.Height());
701 GetMultiplierSolution(ref_constraint_sol);
702 }
703
704 /* Helper routine to reduce code duplication - given a node (which MFEM
705 sometimes calls a "dof"), this returns what normal people call a dof but
706 which MFEM sometimes calls a "vdof" - note that MFEM's naming conventions
707 regarding this are not entirely consistent. In parallel, this always
708 returns the "truedof" in parallel numbering. */
CanonicalNodeNumber(FiniteElementSpace & fespace,int node,bool parallel,int d=0)709 int CanonicalNodeNumber(FiniteElementSpace& fespace,
710 int node, bool parallel, int d=0)
711 {
712 #ifdef MFEM_USE_MPI
713 if (parallel)
714 {
715 ParFiniteElementSpace* pfespace =
716 dynamic_cast<ParFiniteElementSpace*>(&fespace);
717 if (pfespace)
718 {
719 const int vdof = pfespace->DofToVDof(node, d);
720 return pfespace->GetLocalTDofNumber(vdof);
721 }
722 else
723 {
724 MFEM_ABORT("Asked for parallel form of serial object!");
725 return -1;
726 }
727 }
728 else
729 #endif
730 {
731 return fespace.DofToVDof(node, d);
732 }
733 }
734
BuildNormalConstraints(FiniteElementSpace & fespace,Array<int> & constrained_att,Array<int> & constraint_rowstarts,bool parallel)735 SparseMatrix * BuildNormalConstraints(FiniteElementSpace& fespace,
736 Array<int>& constrained_att,
737 Array<int>& constraint_rowstarts,
738 bool parallel)
739 {
740 int dim = fespace.GetVDim();
741
742 // dof_constraint maps a dof (column of the constraint matrix) to
743 // a block-constraint
744 // the indexing is by tdof, but a single tdof uniquely identifies a node
745 // so we only store one tdof independent of dimension
746 std::map<int, int> dof_bconstraint;
747 // constraints[j] is a map from attribute to row number,
748 // the j itself is the index of a block-constraint
749 std::vector<std::map<int, int> > constraints;
750 int n_bconstraints = 0;
751 int n_rows = 0;
752 for (int att : constrained_att)
753 {
754 // identify tdofs on constrained boundary
755 std::set<int> constrained_tdofs;
756 for (int i = 0; i < fespace.GetNBE(); ++i)
757 {
758 if (fespace.GetBdrAttribute(i) == att)
759 {
760 Array<int> nodes;
761 // get nodes on boundary (MFEM sometimes calls these dofs, what
762 // we call dofs it calls vdofs)
763 fespace.GetBdrElementDofs(i, nodes);
764 for (auto k : nodes)
765 {
766 // get the (local) dof number corresponding to
767 // the x-coordinate dof for node k
768 int tdof = CanonicalNodeNumber(fespace, k, parallel);
769 if (tdof >= 0) { constrained_tdofs.insert(tdof); }
770 }
771 }
772 }
773 // fill in the maps identifying which constraints (rows) correspond to
774 // which tdofs
775 for (auto k : constrained_tdofs)
776 {
777 auto it = dof_bconstraint.find(k);
778 if (it == dof_bconstraint.end())
779 {
780 // build new block constraint
781 dof_bconstraint[k] = n_bconstraints++;
782 constraints.emplace_back();
783 constraints.back()[att] = n_rows++;
784 }
785 else
786 {
787 // add tdof to existing block constraint
788 constraints[it->second][att] = n_rows++;
789 }
790 }
791 }
792
793 // reorder so block-constraints eliminated together are grouped together in
794 // adjacent rows
795 {
796 std::map<int, int> reorder_rows;
797 int new_row = 0;
798 constraint_rowstarts.DeleteAll();
799 constraint_rowstarts.Append(0);
800 for (auto& it : dof_bconstraint)
801 {
802 int bconstraint_index = it.second;
803 bool nconstraint = false;
804 for (auto& att_it : constraints[bconstraint_index])
805 {
806 auto rrit = reorder_rows.find(att_it.second);
807 if (rrit == reorder_rows.end())
808 {
809 nconstraint = true;
810 reorder_rows[att_it.second] = new_row++;
811 }
812 }
813 if (nconstraint) { constraint_rowstarts.Append(new_row); }
814 }
815 MFEM_VERIFY(new_row == n_rows, "Remapping failed!");
816 for (auto& constraint_map : constraints)
817 {
818 for (auto& it : constraint_map)
819 {
820 it.second = reorder_rows[it.second];
821 }
822 }
823 }
824
825 SparseMatrix * out = new SparseMatrix(n_rows, fespace.GetTrueVSize());
826
827 // fill in constraint matrix with normal vector information
828 Vector nor(dim);
829 // how many times we have seen a node (key is truek)
830 std::map<int, int> node_visits;
831 for (int i = 0; i < fespace.GetNBE(); ++i)
832 {
833 int att = fespace.GetBdrAttribute(i);
834 if (constrained_att.FindSorted(att) != -1)
835 {
836 ElementTransformation * Tr = fespace.GetBdrElementTransformation(i);
837 const FiniteElement * fe = fespace.GetBE(i);
838 const IntegrationRule& nodes = fe->GetNodes();
839
840 Array<int> dofs;
841 fespace.GetBdrElementDofs(i, dofs);
842 MFEM_VERIFY(dofs.Size() == nodes.Size(),
843 "Something wrong in finite element space!");
844
845 for (int j = 0; j < dofs.Size(); ++j)
846 {
847 Tr->SetIntPoint(&nodes[j]);
848 // the normal returned in the next line is scaled by h, which is
849 // probably what we want in most applications
850 CalcOrtho(Tr->Jacobian(), nor);
851
852 int k = dofs[j];
853 int truek = CanonicalNodeNumber(fespace, k, parallel);
854 if (truek >= 0)
855 {
856 auto nv_it = node_visits.find(truek);
857 if (nv_it == node_visits.end())
858 {
859 node_visits[truek] = 1;
860 }
861 else
862 {
863 node_visits[truek]++;
864 }
865 int visits = node_visits[truek];
866 int bconstraint = dof_bconstraint[truek];
867 int row = constraints[bconstraint][att];
868 for (int d = 0; d < dim; ++d)
869 {
870 int inner_truek = CanonicalNodeNumber(fespace, k,
871 parallel, d);
872 if (visits == 1)
873 {
874 out->Add(row, inner_truek, nor[d]);
875 }
876 else
877 {
878 out->SetColPtr(row);
879 const double pv = out->SearchRow(inner_truek);
880 const double scaling = ((double) (visits - 1)) /
881 ((double) visits);
882 // incremental average, based on how many times
883 // this node has been visited
884 out->Set(row, inner_truek,
885 scaling * pv + (1.0 / visits) * nor[d]);
886 }
887
888 }
889 }
890 }
891 }
892 }
893 out->Finalize();
894
895 return out;
896 }
897
898 #ifdef MFEM_USE_MPI
ParBuildNormalConstraints(ParFiniteElementSpace & fespace,Array<int> & constrained_att,Array<int> & constraint_rowstarts)899 SparseMatrix * ParBuildNormalConstraints(ParFiniteElementSpace& fespace,
900 Array<int>& constrained_att,
901 Array<int>& constraint_rowstarts)
902 {
903 return BuildNormalConstraints(fespace, constrained_att,
904 constraint_rowstarts, true);
905 }
906 #endif
907
908 }
909