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 "../config/config.hpp"
13
14 #ifdef MFEM_USE_MPI
15
16 #include "fem.hpp"
17 #include "../general/sort_pairs.hpp"
18
19 namespace mfem
20 {
21
pAllocMat()22 void ParBilinearForm::pAllocMat()
23 {
24 int nbr_size = pfes->GetFaceNbrVSize();
25
26 if (precompute_sparsity == 0 || fes->GetVDim() > 1)
27 {
28 if (keep_nbr_block)
29 {
30 mat = new SparseMatrix(height + nbr_size, width + nbr_size);
31 }
32 else
33 {
34 mat = new SparseMatrix(height, width + nbr_size);
35 }
36 return;
37 }
38
39 // the sparsity pattern is defined from the map: face->element->dof
40 const Table &lelem_ldof = fes->GetElementToDofTable(); // <-- dofs
41 const Table &nelem_ndof = pfes->face_nbr_element_dof; // <-- vdofs
42 Table elem_dof; // element + nbr-element <---> dof
43 if (nbr_size > 0)
44 {
45 // merge lelem_ldof and nelem_ndof into elem_dof
46 int s1 = lelem_ldof.Size(), s2 = nelem_ndof.Size();
47 const int *I1 = lelem_ldof.GetI(), *J1 = lelem_ldof.GetJ();
48 const int *I2 = nelem_ndof.GetI(), *J2 = nelem_ndof.GetJ();
49 const int nnz1 = I1[s1], nnz2 = I2[s2];
50
51 elem_dof.SetDims(s1 + s2, nnz1 + nnz2);
52
53 int *I = elem_dof.GetI(), *J = elem_dof.GetJ();
54 for (int i = 0; i <= s1; i++)
55 {
56 I[i] = I1[i];
57 }
58 for (int j = 0; j < nnz1; j++)
59 {
60 J[j] = J1[j];
61 }
62 for (int i = 0; i <= s2; i++)
63 {
64 I[s1+i] = I2[i] + nnz1;
65 }
66 for (int j = 0; j < nnz2; j++)
67 {
68 J[nnz1+j] = J2[j] + height;
69 }
70 }
71 // dof_elem x elem_face x face_elem x elem_dof (keep_nbr_block = true)
72 // ldof_lelem x lelem_face x face_elem x elem_dof (keep_nbr_block = false)
73 Table dof_dof;
74 {
75 Table face_dof; // face_elem x elem_dof
76 {
77 Table *face_elem = pfes->GetParMesh()->GetFaceToAllElementTable();
78 if (nbr_size > 0)
79 {
80 mfem::Mult(*face_elem, elem_dof, face_dof);
81 }
82 else
83 {
84 mfem::Mult(*face_elem, lelem_ldof, face_dof);
85 }
86 delete face_elem;
87 if (nbr_size > 0)
88 {
89 elem_dof.Clear();
90 }
91 }
92
93 if (keep_nbr_block)
94 {
95 Table dof_face;
96 Transpose(face_dof, dof_face, height + nbr_size);
97 mfem::Mult(dof_face, face_dof, dof_dof);
98 }
99 else
100 {
101 Table ldof_face;
102 {
103 Table face_ldof;
104 Table *face_lelem = fes->GetMesh()->GetFaceToElementTable();
105 mfem::Mult(*face_lelem, lelem_ldof, face_ldof);
106 delete face_lelem;
107 Transpose(face_ldof, ldof_face, height);
108 }
109 mfem::Mult(ldof_face, face_dof, dof_dof);
110 }
111 }
112
113 int *I = dof_dof.GetI();
114 int *J = dof_dof.GetJ();
115 int nrows = dof_dof.Size();
116 double *data = Memory<double>(I[nrows]);
117
118 mat = new SparseMatrix(I, J, data, nrows, height + nbr_size);
119 *mat = 0.0;
120
121 dof_dof.LoseData();
122 }
123
ParallelAssemble(OperatorHandle & A,SparseMatrix * A_local)124 void ParBilinearForm::ParallelAssemble(OperatorHandle &A, SparseMatrix *A_local)
125 {
126 A.Clear();
127
128 if (A_local == NULL) { return; }
129 MFEM_VERIFY(A_local->Finalized(), "the local matrix must be finalized");
130
131 OperatorHandle dA(A.Type()), Ph(A.Type()), hdA;
132
133 if (interior_face_integs.Size() == 0)
134 {
135 // construct a parallel block-diagonal matrix 'A' based on 'a'
136 dA.MakeSquareBlockDiag(pfes->GetComm(), pfes->GlobalVSize(),
137 pfes->GetDofOffsets(), A_local);
138 }
139 else
140 {
141 // handle the case when 'a' contains off-diagonal
142 int lvsize = pfes->GetVSize();
143 const HYPRE_BigInt *face_nbr_glob_ldof = pfes->GetFaceNbrGlobalDofMap();
144 HYPRE_BigInt ldof_offset = pfes->GetMyDofOffset();
145
146 Array<HYPRE_BigInt> glob_J(A_local->NumNonZeroElems());
147 int *J = A_local->GetJ();
148 for (int i = 0; i < glob_J.Size(); i++)
149 {
150 if (J[i] < lvsize)
151 {
152 glob_J[i] = J[i] + ldof_offset;
153 }
154 else
155 {
156 glob_J[i] = face_nbr_glob_ldof[J[i] - lvsize];
157 }
158 }
159
160 // TODO - construct dA directly in the A format
161 hdA.Reset(
162 new HypreParMatrix(pfes->GetComm(), lvsize, pfes->GlobalVSize(),
163 pfes->GlobalVSize(), A_local->GetI(), glob_J,
164 A_local->GetData(), pfes->GetDofOffsets(),
165 pfes->GetDofOffsets()));
166 // - hdA owns the new HypreParMatrix
167 // - the above constructor copies all input arrays
168 glob_J.DeleteAll();
169 dA.ConvertFrom(hdA);
170 }
171
172 // TODO - assemble the Dof_TrueDof_Matrix directly in the required format?
173 Ph.ConvertFrom(pfes->Dof_TrueDof_Matrix());
174 // TODO: When Ph.Type() == Operator::ANY_TYPE we want to use the Operator
175 // returned by pfes->GetProlongationMatrix(), however that Operator is a
176 // const Operator, so we cannot store it in OperatorHandle. We need a const
177 // version of class OperatorHandle, e.g. ConstOperatorHandle.
178
179 A.MakePtAP(dA, Ph);
180 }
181
ParallelAssemble(SparseMatrix * m)182 HypreParMatrix *ParBilinearForm::ParallelAssemble(SparseMatrix *m)
183 {
184 OperatorHandle Mh(Operator::Hypre_ParCSR);
185 ParallelAssemble(Mh, m);
186 Mh.SetOperatorOwner(false);
187 return Mh.As<HypreParMatrix>();
188 }
189
AssembleSharedFaces(int skip_zeros)190 void ParBilinearForm::AssembleSharedFaces(int skip_zeros)
191 {
192 ParMesh *pmesh = pfes->GetParMesh();
193 FaceElementTransformations *T;
194 Array<int> vdofs1, vdofs2, vdofs_all;
195 DenseMatrix elemmat;
196
197 int nfaces = pmesh->GetNSharedFaces();
198 for (int i = 0; i < nfaces; i++)
199 {
200 T = pmesh->GetSharedFaceTransformations(i);
201 int Elem2NbrNo = T->Elem2No - pmesh->GetNE();
202 pfes->GetElementVDofs(T->Elem1No, vdofs1);
203 pfes->GetFaceNbrElementVDofs(Elem2NbrNo, vdofs2);
204 vdofs1.Copy(vdofs_all);
205 for (int j = 0; j < vdofs2.Size(); j++)
206 {
207 if (vdofs2[j] >= 0)
208 {
209 vdofs2[j] += height;
210 }
211 else
212 {
213 vdofs2[j] -= height;
214 }
215 }
216 vdofs_all.Append(vdofs2);
217 for (int k = 0; k < interior_face_integs.Size(); k++)
218 {
219 interior_face_integs[k]->
220 AssembleFaceMatrix(*pfes->GetFE(T->Elem1No),
221 *pfes->GetFaceNbrFE(Elem2NbrNo),
222 *T, elemmat);
223 if (keep_nbr_block)
224 {
225 mat->AddSubMatrix(vdofs_all, vdofs_all, elemmat, skip_zeros);
226 }
227 else
228 {
229 mat->AddSubMatrix(vdofs1, vdofs_all, elemmat, skip_zeros);
230 }
231 }
232 }
233 }
234
Assemble(int skip_zeros)235 void ParBilinearForm::Assemble(int skip_zeros)
236 {
237 if (interior_face_integs.Size())
238 {
239 pfes->ExchangeFaceNbrData();
240 if (!ext && mat == NULL)
241 {
242 pAllocMat();
243 }
244 }
245
246 BilinearForm::Assemble(skip_zeros);
247
248 if (!ext && interior_face_integs.Size() > 0)
249 {
250 AssembleSharedFaces(skip_zeros);
251 }
252 }
253
AssembleDiagonal(Vector & diag) const254 void ParBilinearForm::AssembleDiagonal(Vector &diag) const
255 {
256 MFEM_ASSERT(diag.Size() == fes->GetTrueVSize(),
257 "Vector for holding diagonal has wrong size!");
258 const Operator *P = fes->GetProlongationMatrix();
259 if (!ext)
260 {
261 MFEM_ASSERT(p_mat.Ptr(), "the ParBilinearForm is not assembled!");
262 p_mat->AssembleDiagonal(diag); // TODO: add support for PETSc matrices
263 return;
264 }
265 // Here, we have extension, ext.
266 if (IsIdentityProlongation(P))
267 {
268 ext->AssembleDiagonal(diag);
269 return;
270 }
271 // Here, we have extension, ext, and parallel/conforming prolongation, P.
272 Vector local_diag(P->Height());
273 ext->AssembleDiagonal(local_diag);
274 if (fes->Conforming())
275 {
276 P->MultTranspose(local_diag, diag);
277 return;
278 }
279 // For an AMR mesh, a convergent diagonal is assembled with |P^T| d_l,
280 // where |P^T| has the entry-wise absolute values of the conforming
281 // prolongation transpose operator.
282 const HypreParMatrix *HP = dynamic_cast<const HypreParMatrix*>(P);
283 if (HP)
284 {
285 HP->AbsMultTranspose(1.0, local_diag, 0.0, diag);
286 }
287 else
288 {
289 MFEM_ABORT("unsupported prolongation matrix type.");
290 }
291 }
292
293 void ParBilinearForm
ParallelEliminateEssentialBC(const Array<int> & bdr_attr_is_ess,HypreParMatrix & A,const HypreParVector & X,HypreParVector & B) const294 ::ParallelEliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
295 HypreParMatrix &A, const HypreParVector &X,
296 HypreParVector &B) const
297 {
298 Array<int> dof_list;
299
300 pfes->GetEssentialTrueDofs(bdr_attr_is_ess, dof_list);
301
302 // do the parallel elimination
303 A.EliminateRowsCols(dof_list, X, B);
304 }
305
306 HypreParMatrix *ParBilinearForm::
ParallelEliminateEssentialBC(const Array<int> & bdr_attr_is_ess,HypreParMatrix & A) const307 ParallelEliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
308 HypreParMatrix &A) const
309 {
310 Array<int> dof_list;
311
312 pfes->GetEssentialTrueDofs(bdr_attr_is_ess, dof_list);
313
314 return A.EliminateRowsCols(dof_list);
315 }
316
TrueAddMult(const Vector & x,Vector & y,const double a) const317 void ParBilinearForm::TrueAddMult(const Vector &x, Vector &y, const double a)
318 const
319 {
320 MFEM_VERIFY(interior_face_integs.Size() == 0,
321 "the case of interior face integrators is not"
322 " implemented");
323
324 if (X.ParFESpace() != pfes)
325 {
326 X.SetSpace(pfes);
327 Y.SetSpace(pfes);
328 }
329
330 X.Distribute(&x);
331 if (ext)
332 {
333 ext->Mult(X, Y);
334 }
335 else
336 {
337 mat->Mult(X, Y);
338 }
339 pfes->Dof_TrueDof_Matrix()->MultTranspose(a, Y, 1.0, y);
340 }
341
FormLinearSystem(const Array<int> & ess_tdof_list,Vector & x,Vector & b,OperatorHandle & A,Vector & X,Vector & B,int copy_interior)342 void ParBilinearForm::FormLinearSystem(
343 const Array<int> &ess_tdof_list, Vector &x, Vector &b,
344 OperatorHandle &A, Vector &X, Vector &B, int copy_interior)
345 {
346 if (ext)
347 {
348 ext->FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
349 return;
350 }
351
352 // Finish the matrix assembly and perform BC elimination, storing the
353 // eliminated part of the matrix.
354 FormSystemMatrix(ess_tdof_list, A);
355
356 const Operator &P = *pfes->GetProlongationMatrix();
357 const SparseMatrix &R = *pfes->GetRestrictionMatrix();
358
359 // Transform the system and perform the elimination in B, based on the
360 // essential BC values from x. Restrict the BC part of x in X, and set the
361 // non-BC part to zero. Since there is no good initial guess for the Lagrange
362 // multipliers, set X = 0.0 for hybridization.
363 if (static_cond)
364 {
365 // Schur complement reduction to the exposed dofs
366 static_cond->ReduceSystem(x, b, X, B, copy_interior);
367 }
368 else if (hybridization)
369 {
370 // Reduction to the Lagrange multipliers system
371 HypreParVector true_X(pfes), true_B(pfes);
372 P.MultTranspose(b, true_B);
373 R.Mult(x, true_X);
374 p_mat.EliminateBC(p_mat_e, ess_tdof_list, true_X, true_B);
375 R.MultTranspose(true_B, b);
376 hybridization->ReduceRHS(true_B, B);
377 X.SetSize(B.Size());
378 X = 0.0;
379 }
380 else
381 {
382 // Variational restriction with P
383 X.SetSize(P.Width());
384 B.SetSize(X.Size());
385 P.MultTranspose(b, B);
386 R.Mult(x, X);
387 p_mat.EliminateBC(p_mat_e, ess_tdof_list, X, B);
388 if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
389 }
390 }
391
EliminateVDofsInRHS(const Array<int> & vdofs,const Vector & x,Vector & b)392 void ParBilinearForm::EliminateVDofsInRHS(
393 const Array<int> &vdofs, const Vector &x, Vector &b)
394 {
395 p_mat.EliminateBC(p_mat_e, vdofs, x, b);
396 }
397
FormSystemMatrix(const Array<int> & ess_tdof_list,OperatorHandle & A)398 void ParBilinearForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
399 OperatorHandle &A)
400 {
401 if (ext)
402 {
403 ext->FormSystemMatrix(ess_tdof_list, A);
404 return;
405 }
406
407 // Finish the matrix assembly and perform BC elimination, storing the
408 // eliminated part of the matrix.
409 if (static_cond)
410 {
411 if (!static_cond->HasEliminatedBC())
412 {
413 static_cond->SetEssentialTrueDofs(ess_tdof_list);
414 static_cond->Finalize();
415 static_cond->EliminateReducedTrueDofs(Matrix::DIAG_ONE);
416 }
417 static_cond->GetParallelMatrix(A);
418 }
419 else
420 {
421 if (mat)
422 {
423 const int remove_zeros = 0;
424 Finalize(remove_zeros);
425 MFEM_VERIFY(p_mat.Ptr() == NULL && p_mat_e.Ptr() == NULL,
426 "The ParBilinearForm must be updated with Update() before "
427 "re-assembling the ParBilinearForm.");
428 ParallelAssemble(p_mat, mat);
429 delete mat;
430 mat = NULL;
431 delete mat_e;
432 mat_e = NULL;
433 p_mat_e.EliminateRowsCols(p_mat, ess_tdof_list);
434 }
435 if (hybridization)
436 {
437 hybridization->GetParallelMatrix(A);
438 }
439 else
440 {
441 A = p_mat;
442 }
443 }
444 }
445
RecoverFEMSolution(const Vector & X,const Vector & b,Vector & x)446 void ParBilinearForm::RecoverFEMSolution(
447 const Vector &X, const Vector &b, Vector &x)
448 {
449 if (ext)
450 {
451 ext->RecoverFEMSolution(X, b, x);
452 return;
453 }
454
455 const Operator &P = *pfes->GetProlongationMatrix();
456
457 if (static_cond)
458 {
459 // Private dofs back solve
460 static_cond->ComputeSolution(b, X, x);
461 }
462 else if (hybridization)
463 {
464 // Primal unknowns recovery
465 HypreParVector true_X(pfes), true_B(pfes);
466 P.MultTranspose(b, true_B);
467 const SparseMatrix &R = *pfes->GetRestrictionMatrix();
468 R.Mult(x, true_X); // get essential b.c. from x
469 hybridization->ComputeSolution(true_B, X, true_X);
470 x.SetSize(P.Height());
471 P.Mult(true_X, x);
472 }
473 else
474 {
475 // Apply conforming prolongation
476 x.SetSize(P.Height(), GetHypreMemoryType());
477 P.Mult(X, x);
478 }
479 }
480
Update(FiniteElementSpace * nfes)481 void ParBilinearForm::Update(FiniteElementSpace *nfes)
482 {
483 BilinearForm::Update(nfes);
484
485 if (nfes)
486 {
487 pfes = dynamic_cast<ParFiniteElementSpace *>(nfes);
488 MFEM_VERIFY(pfes != NULL, "nfes must be a ParFiniteElementSpace!");
489 }
490
491 p_mat.Clear();
492 p_mat_e.Clear();
493 }
494
495
ParallelAssemble()496 HypreParMatrix *ParMixedBilinearForm::ParallelAssemble()
497 {
498 // construct the block-diagonal matrix A
499 HypreParMatrix *A =
500 new HypreParMatrix(trial_pfes->GetComm(),
501 test_pfes->GlobalVSize(),
502 trial_pfes->GlobalVSize(),
503 test_pfes->GetDofOffsets(),
504 trial_pfes->GetDofOffsets(),
505 mat);
506
507 HypreParMatrix *rap = RAP(test_pfes->Dof_TrueDof_Matrix(), A,
508 trial_pfes->Dof_TrueDof_Matrix());
509
510 delete A;
511
512 return rap;
513 }
514
ParallelAssemble(OperatorHandle & A)515 void ParMixedBilinearForm::ParallelAssemble(OperatorHandle &A)
516 {
517 // construct the rectangular block-diagonal matrix dA
518 OperatorHandle dA(A.Type());
519 dA.MakeRectangularBlockDiag(trial_pfes->GetComm(),
520 test_pfes->GlobalVSize(),
521 trial_pfes->GlobalVSize(),
522 test_pfes->GetDofOffsets(),
523 trial_pfes->GetDofOffsets(),
524 mat);
525
526 OperatorHandle P_test(A.Type()), P_trial(A.Type());
527
528 // TODO - construct the Dof_TrueDof_Matrix directly in the required format.
529 P_test.ConvertFrom(test_pfes->Dof_TrueDof_Matrix());
530 P_trial.ConvertFrom(trial_pfes->Dof_TrueDof_Matrix());
531
532 A.MakeRAP(P_test, dA, P_trial);
533 }
534
535 /// Compute y += a (P^t A P) x, where x and y are vectors on the true dofs
TrueAddMult(const Vector & x,Vector & y,const double a) const536 void ParMixedBilinearForm::TrueAddMult(const Vector &x, Vector &y,
537 const double a) const
538 {
539 if (X.ParFESpace() != trial_pfes)
540 {
541 X.SetSpace(trial_pfes);
542 Y.SetSpace(test_pfes);
543 }
544
545 X.Distribute(&x);
546 mat->Mult(X, Y);
547 test_pfes->Dof_TrueDof_Matrix()->MultTranspose(a, Y, 1.0, y);
548 }
549
FormRectangularSystemMatrix(const Array<int> & trial_tdof_list,const Array<int> & test_tdof_list,OperatorHandle & A)550 void ParMixedBilinearForm::FormRectangularSystemMatrix(
551 const Array<int>
552 &trial_tdof_list,
553 const Array<int> &test_tdof_list,
554 OperatorHandle &A)
555 {
556 if (ext)
557 {
558 ext->FormRectangularSystemOperator(trial_tdof_list, test_tdof_list, A);
559 return;
560 }
561
562 if (mat)
563 {
564 Finalize();
565 ParallelAssemble(p_mat);
566 delete mat;
567 mat = NULL;
568 delete mat_e;
569 mat_e = NULL;
570 HypreParMatrix *temp =
571 p_mat.As<HypreParMatrix>()->EliminateCols(trial_tdof_list);
572 p_mat.As<HypreParMatrix>()->EliminateRows(test_tdof_list);
573 p_mat_e.Reset(temp, true);
574 }
575
576 A = p_mat;
577 }
578
FormRectangularLinearSystem(const Array<int> & trial_tdof_list,const Array<int> & test_tdof_list,Vector & x,Vector & b,OperatorHandle & A,Vector & X,Vector & B)579 void ParMixedBilinearForm::FormRectangularLinearSystem(
580 const Array<int>
581 &trial_tdof_list,
582 const Array<int> &test_tdof_list, Vector &x,
583 Vector &b, OperatorHandle &A, Vector &X,
584 Vector &B)
585 {
586 if (ext)
587 {
588 ext->FormRectangularLinearSystem(trial_tdof_list, test_tdof_list,
589 x, b, A, X, B);
590 return;
591 }
592
593 FormRectangularSystemMatrix(trial_tdof_list, test_tdof_list, A);
594
595 const Operator *test_P = test_pfes->GetProlongationMatrix();
596 const SparseMatrix *trial_R = trial_pfes->GetRestrictionMatrix();
597
598 X.SetSize(trial_pfes->TrueVSize());
599 B.SetSize(test_pfes->TrueVSize());
600 test_P->MultTranspose(b, B);
601 trial_R->Mult(x, X);
602
603 p_mat_e.As<HypreParMatrix>()->Mult(-1.0, X, 1.0, B);
604 B.SetSubVector(test_tdof_list, 0.0);
605 }
606
ParallelAssemble() const607 HypreParMatrix* ParDiscreteLinearOperator::ParallelAssemble() const
608 {
609 MFEM_ASSERT(mat, "Matrix is not assembled");
610 MFEM_ASSERT(mat->Finalized(), "Matrix is not finalized");
611 SparseMatrix* RA = mfem::Mult(*range_fes->GetRestrictionMatrix(), *mat);
612 HypreParMatrix* P = domain_fes->Dof_TrueDof_Matrix();
613 HypreParMatrix* RAP = P->LeftDiagMult(*RA, range_fes->GetTrueDofOffsets());
614 delete RA;
615 return RAP;
616 }
617
ParallelAssemble(OperatorHandle & A)618 void ParDiscreteLinearOperator::ParallelAssemble(OperatorHandle &A)
619 {
620 // construct the rectangular block-diagonal matrix dA
621 OperatorHandle dA(A.Type());
622 dA.MakeRectangularBlockDiag(domain_fes->GetComm(),
623 range_fes->GlobalVSize(),
624 domain_fes->GlobalVSize(),
625 range_fes->GetDofOffsets(),
626 domain_fes->GetDofOffsets(),
627 mat);
628
629 OperatorHandle R_test_transpose(A.Type()), P_trial(A.Type());
630
631 // TODO - construct the Dof_TrueDof_Matrix directly in the required format.
632 R_test_transpose.ConvertFrom(range_fes->Dof_TrueDof_Matrix());
633 P_trial.ConvertFrom(domain_fes->Dof_TrueDof_Matrix());
634
635 A.MakeRAP(R_test_transpose, dA, P_trial);
636 }
637
FormRectangularSystemMatrix(OperatorHandle & A)638 void ParDiscreteLinearOperator::FormRectangularSystemMatrix(OperatorHandle &A)
639 {
640 if (ext)
641 {
642 Array<int> empty;
643 ext->FormRectangularSystemOperator(empty, empty, A);
644 return;
645 }
646
647 mfem_error("not implemented!");
648 }
649
GetParBlocks(Array2D<HypreParMatrix * > & blocks) const650 void ParDiscreteLinearOperator::GetParBlocks(Array2D<HypreParMatrix *> &blocks)
651 const
652 {
653 MFEM_VERIFY(mat->Finalized(), "Local matrix needs to be finalized for "
654 "GetParBlocks");
655
656 HypreParMatrix* RLP = ParallelAssemble();
657
658 blocks.SetSize(range_fes->GetVDim(), domain_fes->GetVDim());
659
660 RLP->GetBlocks(blocks,
661 range_fes->GetOrdering() == Ordering::byVDIM,
662 domain_fes->GetOrdering() == Ordering::byVDIM);
663
664 delete RLP;
665 }
666
667 }
668
669 #endif
670