MFEM  v4.6.0
Finite element discretization library
pbilinearform.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 "../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 
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 
125  bool steal_loc_A)
126 {
127  ParFiniteElementSpace &pfespace = *ParFESpace();
128 
129  // Create a block diagonal parallel matrix
131  A_diag.MakeSquareBlockDiag(pfespace.GetComm(),
132  pfespace.GlobalVSize(),
133  pfespace.GetDofOffsets(),
134  &loc_A);
135 
136  // Parallel matrix assembly using P^t A P (if needed)
138  {
139  A_diag.SetOperatorOwner(false);
140  A.Reset(A_diag.As<HypreParMatrix>());
141  if (steal_loc_A)
142  {
143  HypreStealOwnership(*A.As<HypreParMatrix>(), loc_A);
144  }
145  }
146  else
147  {
149  P.ConvertFrom(pfespace.Dof_TrueDof_Matrix());
150  A.MakePtAP(A_diag, P);
151  }
152 }
153 
155 {
156  A.Clear();
157 
158  if (A_local == NULL) { return; }
159  MFEM_VERIFY(A_local->Finalized(), "the local matrix must be finalized");
160 
161  OperatorHandle dA(A.Type()), Ph(A.Type()), hdA;
162 
163  if (interior_face_integs.Size() == 0)
164  {
165  // construct a parallel block-diagonal matrix 'A' based on 'a'
167  pfes->GetDofOffsets(), A_local);
168  }
169  else
170  {
171  // handle the case when 'a' contains off-diagonal
172  int lvsize = pfes->GetVSize();
173  const HYPRE_BigInt *face_nbr_glob_ldof = pfes->GetFaceNbrGlobalDofMap();
174  HYPRE_BigInt ldof_offset = pfes->GetMyDofOffset();
175 
176  Array<HYPRE_BigInt> glob_J(A_local->NumNonZeroElems());
177  int *J = A_local->GetJ();
178  for (int i = 0; i < glob_J.Size(); i++)
179  {
180  if (J[i] < lvsize)
181  {
182  glob_J[i] = J[i] + ldof_offset;
183  }
184  else
185  {
186  glob_J[i] = face_nbr_glob_ldof[J[i] - lvsize];
187  }
188  }
189 
190  // TODO - construct dA directly in the A format
191  hdA.Reset(
192  new HypreParMatrix(pfes->GetComm(), lvsize, pfes->GlobalVSize(),
193  pfes->GlobalVSize(), A_local->GetI(), glob_J,
194  A_local->GetData(), pfes->GetDofOffsets(),
195  pfes->GetDofOffsets()));
196  // - hdA owns the new HypreParMatrix
197  // - the above constructor copies all input arrays
198  glob_J.DeleteAll();
199  dA.ConvertFrom(hdA);
200  }
201 
202  // TODO - assemble the Dof_TrueDof_Matrix directly in the required format?
203  Ph.ConvertFrom(pfes->Dof_TrueDof_Matrix());
204  // TODO: When Ph.Type() == Operator::ANY_TYPE we want to use the Operator
205  // returned by pfes->GetProlongationMatrix(), however that Operator is a
206  // const Operator, so we cannot store it in OperatorHandle. We need a const
207  // version of class OperatorHandle, e.g. ConstOperatorHandle.
208 
209  A.MakePtAP(dA, Ph);
210 }
211 
213 {
215  ParallelAssemble(Mh, m);
216  Mh.SetOperatorOwner(false);
217  return Mh.As<HypreParMatrix>();
218 }
219 
221 {
222  ParMesh *pmesh = pfes->GetParMesh();
224  Array<int> vdofs1, vdofs2, vdofs_all;
226 
227  int nfaces = pmesh->GetNSharedFaces();
228  for (int i = 0; i < nfaces; i++)
229  {
230  T = pmesh->GetSharedFaceTransformations(i);
231  int Elem2NbrNo = T->Elem2No - pmesh->GetNE();
232  pfes->GetElementVDofs(T->Elem1No, vdofs1);
233  pfes->GetFaceNbrElementVDofs(Elem2NbrNo, vdofs2);
234  vdofs1.Copy(vdofs_all);
235  for (int j = 0; j < vdofs2.Size(); j++)
236  {
237  if (vdofs2[j] >= 0)
238  {
239  vdofs2[j] += height;
240  }
241  else
242  {
243  vdofs2[j] -= height;
244  }
245  }
246  vdofs_all.Append(vdofs2);
247  for (int k = 0; k < interior_face_integs.Size(); k++)
248  {
250  AssembleFaceMatrix(*pfes->GetFE(T->Elem1No),
251  *pfes->GetFaceNbrFE(Elem2NbrNo),
252  *T, elemmat);
253  if (keep_nbr_block)
254  {
255  mat->AddSubMatrix(vdofs_all, vdofs_all, elemmat, skip_zeros);
256  }
257  else
258  {
259  mat->AddSubMatrix(vdofs1, vdofs_all, elemmat, skip_zeros);
260  }
261  }
262  }
263 }
264 
265 void ParBilinearForm::Assemble(int skip_zeros)
266 {
267  if (interior_face_integs.Size())
268  {
270  if (!ext && mat == NULL)
271  {
272  pAllocMat();
273  }
274  }
275 
276  BilinearForm::Assemble(skip_zeros);
277 
278  if (!ext && interior_face_integs.Size() > 0)
279  {
280  AssembleSharedFaces(skip_zeros);
281  }
282 }
283 
285 {
286  MFEM_ASSERT(diag.Size() == fes->GetTrueVSize(),
287  "Vector for holding diagonal has wrong size!");
288  const Operator *P = fes->GetProlongationMatrix();
289  if (!ext)
290  {
291  MFEM_ASSERT(p_mat.Ptr(), "the ParBilinearForm is not assembled!");
292  p_mat->AssembleDiagonal(diag); // TODO: add support for PETSc matrices
293  return;
294  }
295  // Here, we have extension, ext.
296  if (IsIdentityProlongation(P))
297  {
298  ext->AssembleDiagonal(diag);
299  return;
300  }
301  // Here, we have extension, ext, and parallel/conforming prolongation, P.
302  Vector local_diag(P->Height());
303  ext->AssembleDiagonal(local_diag);
304  if (fes->Conforming())
305  {
306  P->MultTranspose(local_diag, diag);
307  return;
308  }
309  // For an AMR mesh, a convergent diagonal is assembled with |P^T| d_l,
310  // where |P^T| has the entry-wise absolute values of the conforming
311  // prolongation transpose operator.
312  const HypreParMatrix *HP = dynamic_cast<const HypreParMatrix*>(P);
313  if (HP)
314  {
315  HP->AbsMultTranspose(1.0, local_diag, 0.0, diag);
316  }
317  else
318  {
319  MFEM_ABORT("unsupported prolongation matrix type.");
320  }
321 }
322 
323 void ParBilinearForm
325  HypreParMatrix &A, const HypreParVector &X,
326  HypreParVector &B) const
327 {
328  Array<int> dof_list;
329 
330  pfes->GetEssentialTrueDofs(bdr_attr_is_ess, dof_list);
331 
332  // do the parallel elimination
333  A.EliminateRowsCols(dof_list, X, B);
334 }
335 
338  HypreParMatrix &A) const
339 {
340  Array<int> dof_list;
341 
342  pfes->GetEssentialTrueDofs(bdr_attr_is_ess, dof_list);
343 
344  return A.EliminateRowsCols(dof_list);
345 }
346 
347 void ParBilinearForm::TrueAddMult(const Vector &x, Vector &y, const double a)
348 const
349 {
350  const Operator *P = pfes->GetProlongationMatrix();
351  Xaux.SetSize(P->Height());
352  Yaux.SetSize(P->Height());
353  Ytmp.SetSize(P->Width());
354 
355  P->Mult(x, Xaux);
356  if (ext)
357  {
358  ext->Mult(Xaux, Yaux);
359  }
360  else
361  {
362  MFEM_VERIFY(interior_face_integs.Size() == 0,
363  "the case of interior face integrators is not"
364  " implemented");
365  mat->Mult(Xaux, Yaux);
366  }
367  P->MultTranspose(Yaux, Ytmp);
368  y.Add(a, Ytmp);
369 }
370 
372  const Array<int> &ess_tdof_list, Vector &x, Vector &b,
373  OperatorHandle &A, Vector &X, Vector &B, int copy_interior)
374 {
375  if (ext)
376  {
377  ext->FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
378  return;
379  }
380 
381  // Finish the matrix assembly and perform BC elimination, storing the
382  // eliminated part of the matrix.
383  FormSystemMatrix(ess_tdof_list, A);
384 
385  const Operator &P = *pfes->GetProlongationMatrix();
386  const SparseMatrix &R = *pfes->GetRestrictionMatrix();
387 
388  // Transform the system and perform the elimination in B, based on the
389  // essential BC values from x. Restrict the BC part of x in X, and set the
390  // non-BC part to zero. Since there is no good initial guess for the Lagrange
391  // multipliers, set X = 0.0 for hybridization.
392  if (static_cond)
393  {
394  // Schur complement reduction to the exposed dofs
395  static_cond->ReduceSystem(x, b, X, B, copy_interior);
396  }
397  else if (hybridization)
398  {
399  // Reduction to the Lagrange multipliers system
400  HypreParVector true_X(pfes), true_B(pfes);
401  P.MultTranspose(b, true_B);
402  R.Mult(x, true_X);
403  p_mat.EliminateBC(p_mat_e, ess_tdof_list, true_X, true_B);
404  R.MultTranspose(true_B, b);
405  hybridization->ReduceRHS(true_B, B);
406  X.SetSize(B.Size());
407  X = 0.0;
408  }
409  else
410  {
411  // Variational restriction with P
412  X.SetSize(P.Width());
413  B.SetSize(X.Size());
414  P.MultTranspose(b, B);
415  R.Mult(x, X);
416  p_mat.EliminateBC(p_mat_e, ess_tdof_list, X, B);
417  if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
418  }
419 }
420 
422  const Array<int> &vdofs, const Vector &x, Vector &b)
423 {
425 }
426 
428  OperatorHandle &A)
429 {
430  if (ext)
431  {
432  ext->FormSystemMatrix(ess_tdof_list, A);
433  return;
434  }
435 
436  // Finish the matrix assembly and perform BC elimination, storing the
437  // eliminated part of the matrix.
438  if (static_cond)
439  {
441  {
442  static_cond->SetEssentialTrueDofs(ess_tdof_list);
445  }
447  }
448  else
449  {
450  if (mat)
451  {
452  const int remove_zeros = 0;
453  Finalize(remove_zeros);
454  MFEM_VERIFY(p_mat.Ptr() == NULL && p_mat_e.Ptr() == NULL,
455  "The ParBilinearForm must be updated with Update() before "
456  "re-assembling the ParBilinearForm.");
458  delete mat;
459  mat = NULL;
460  delete mat_e;
461  mat_e = NULL;
462  p_mat_e.EliminateRowsCols(p_mat, ess_tdof_list);
463  }
464  if (hybridization)
465  {
467  }
468  else
469  {
470  A = p_mat;
471  }
472  }
473 }
474 
476  const Vector &X, const Vector &b, Vector &x)
477 {
478  if (ext)
479  {
480  ext->RecoverFEMSolution(X, b, x);
481  return;
482  }
483 
484  const Operator &P = *pfes->GetProlongationMatrix();
485 
486  if (static_cond)
487  {
488  // Private dofs back solve
489  static_cond->ComputeSolution(b, X, x);
490  }
491  else if (hybridization)
492  {
493  // Primal unknowns recovery
494  HypreParVector true_X(pfes), true_B(pfes);
495  P.MultTranspose(b, true_B);
496  const SparseMatrix &R = *pfes->GetRestrictionMatrix();
497  R.Mult(x, true_X); // get essential b.c. from x
498  hybridization->ComputeSolution(true_B, X, true_X);
499  x.SetSize(P.Height());
500  P.Mult(true_X, x);
501  }
502  else
503  {
504  // Apply conforming prolongation
506  P.Mult(X, x);
507  }
508 }
509 
511 {
512  BilinearForm::Update(nfes);
513 
514  if (nfes)
515  {
516  pfes = dynamic_cast<ParFiniteElementSpace *>(nfes);
517  MFEM_VERIFY(pfes != NULL, "nfes must be a ParFiniteElementSpace!");
518  }
519 
520  p_mat.Clear();
521  p_mat_e.Clear();
522 }
523 
524 
526 {
527  // construct the block-diagonal matrix A
528  HypreParMatrix *A =
534  mat);
535 
538 
539  delete A;
540 
541  return rap;
542 }
543 
545 {
546  // construct the rectangular block-diagonal matrix dA
547  OperatorHandle dA(A.Type());
553  mat);
554 
555  OperatorHandle P_test(A.Type()), P_trial(A.Type());
556 
557  // TODO - construct the Dof_TrueDof_Matrix directly in the required format.
559  P_trial.ConvertFrom(trial_pfes->Dof_TrueDof_Matrix());
560 
561  A.MakeRAP(P_test, dA, P_trial);
562 }
563 
564 /// Compute y += a (P^t A P) x, where x and y are vectors on the true dofs
566  const double a) const
567 {
568  if (Xaux.ParFESpace() != trial_pfes)
569  {
572  }
573 
574  Xaux.Distribute(&x);
575  mat->Mult(Xaux, Yaux);
577 }
578 
580  const Array<int>
581  &trial_tdof_list,
582  const Array<int> &test_tdof_list,
583  OperatorHandle &A)
584 {
585  if (ext)
586  {
587  ext->FormRectangularSystemOperator(trial_tdof_list, test_tdof_list, A);
588  return;
589  }
590 
591  if (mat)
592  {
593  Finalize();
595  delete mat;
596  mat = NULL;
597  delete mat_e;
598  mat_e = NULL;
599  HypreParMatrix *temp =
600  p_mat.As<HypreParMatrix>()->EliminateCols(trial_tdof_list);
601  p_mat.As<HypreParMatrix>()->EliminateRows(test_tdof_list);
602  p_mat_e.Reset(temp, true);
603  }
604 
605  A = p_mat;
606 }
607 
609  const Array<int>
610  &trial_tdof_list,
611  const Array<int> &test_tdof_list, Vector &x,
612  Vector &b, OperatorHandle &A, Vector &X,
613  Vector &B)
614 {
615  if (ext)
616  {
617  ext->FormRectangularLinearSystem(trial_tdof_list, test_tdof_list,
618  x, b, A, X, B);
619  return;
620  }
621 
622  FormRectangularSystemMatrix(trial_tdof_list, test_tdof_list, A);
623 
624  const Operator *test_P = test_pfes->GetProlongationMatrix();
625  const SparseMatrix *trial_R = trial_pfes->GetRestrictionMatrix();
626 
629  test_P->MultTranspose(b, B);
630  trial_R->Mult(x, X);
631 
632  p_mat_e.As<HypreParMatrix>()->Mult(-1.0, X, 1.0, B);
633  B.SetSubVector(test_tdof_list, 0.0);
634 }
635 
637 {
638  MFEM_ASSERT(mat, "Matrix is not assembled");
639  MFEM_ASSERT(mat->Finalized(), "Matrix is not finalized");
642  HypreParMatrix* RAP = P->LeftDiagMult(*RA, range_fes->GetTrueDofOffsets());
643  delete RA;
644  return RAP;
645 }
646 
648 {
649  // construct the rectangular block-diagonal matrix dA
650  OperatorHandle dA(A.Type());
656  mat);
657 
659  OperatorHandle R_test_transpose(A.Type());
660  R_test_transpose.MakeRectangularBlockDiag(range_fes->GetComm(),
665  Rt);
666 
667  // TODO - construct the Dof_TrueDof_Matrix directly in the required format.
668  OperatorHandle P_trial(A.Type());
670 
671  A.MakeRAP(R_test_transpose, dA, P_trial);
672  delete Rt;
673 }
674 
676 {
677  if (ext)
678  {
679  Array<int> empty;
680  ext->FormRectangularSystemOperator(empty, empty, A);
681  return;
682  }
683 
684  mfem_error("not implemented!");
685 }
686 
688 const
689 {
690  MFEM_VERIFY(mat->Finalized(), "Local matrix needs to be finalized for "
691  "GetParBlocks");
692 
694 
695  blocks.SetSize(range_fes->GetVDim(), domain_fes->GetVDim());
696 
697  RLP->GetBlocks(blocks,
700 
701  delete RLP;
702 }
703 
704 }
705 
706 #endif
MemoryType GetHypreMemoryType()
The MemoryType used by MFEM when allocating arrays for Hypre objects.
Definition: hypre.hpp:149
virtual void FormRectangularSystemMatrix(OperatorHandle &A)
Return in A a parallel (on truedofs) version of this operator.
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
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Definition: sparsemat.cpp:733
void EliminateRowsCols(const Array< int > &rows_cols, const HypreParVector &X, HypreParVector &B)
Definition: hypre.cpp:2276
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:605
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Definition: pfespace.cpp:1031
virtual int NumNonZeroElems() const
Returns the number of the nonzero elements in the matrix.
Definition: sparsemat.cpp:1588
void ConvertFrom(OperatorHandle &A)
Convert the given OperatorHandle A to the currently set type id.
Definition: handle.cpp:200
ParFiniteElementSpace * ParFESpace() const
Return the parallel FE space associated with the ParBilinearForm.
int * GetJ()
Definition: table.hpp:114
virtual const Operator * GetProlongationMatrix() const
The returned Operator is owned by the FiniteElementSpace.
Definition: pfespace.cpp:1152
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
void ComputeSolution(const Vector &b, const Vector &sc_sol, Vector &sol) const
Definition: staticcond.cpp:476
virtual void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
Form the linear system A X = B, corresponding to this bilinear form and the linear form b(...
const FiniteElement * GetFaceNbrFE(int i) const
Definition: pfespace.cpp:1510
int TrueVSize() const
Obsolete, kept for backward compatibility.
Definition: pfespace.hpp:424
BilinearFormExtension * ext
Extension for supporting Full Assembly (FA), Element Assembly (EA), Partial Assembly (PA)...
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:517
void TrueAddMult(const Vector &x, Vector &y, const double a=1.0) const
Compute y += a (P^t A P) x, where x and y are vectors on the true dofs.
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
HypreParMatrix & GetParallelMatrix()
Return the parallel hybridized matrix.
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
int * GetJ()
Return the array J.
Definition: sparsemat.hpp:227
void SetDims(int rows, int nnz)
Definition: table.cpp:140
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
void ParallelRAP(SparseMatrix &loc_A, OperatorHandle &A, bool steal_loc_A=false)
Compute parallel RAP operator and store it in A as a HypreParMatrix.
HypreParMatrix * ParallelAssemble() const
Returns the matrix "assembled" on the true dofs.
A specialized ElementTransformation class representing a face and its two neighboring elements...
Definition: eltrans.hpp:480
int Size() const
Returns the size of the vector.
Definition: vector.hpp:197
Array< BilinearFormIntegrator * > interior_face_integs
Set of interior face Integrators to be applied.
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
bool HasEliminatedBC() const
Return true if essential boundary conditions have been eliminated from the Schur complement matrix...
Definition: staticcond.hpp:142
int * GetI()
Return the array I.
Definition: sparsemat.hpp:222
HYPRE_BigInt GetMyDofOffset() const
Definition: pfespace.cpp:1142
void EliminateReducedTrueDofs(const Array< int > &ess_rtdof_list, Matrix::DiagonalPolicy dpolicy)
Eliminate the given reduced true dofs from the Schur complement matrix S.
Definition: staticcond.cpp:286
OperatorHandle p_mat_e
void GetParBlocks(Array2D< HypreParMatrix *> &blocks) const
void LoseData()
Call this if data has been stolen.
Definition: table.hpp:180
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
virtual void SetSpace(FiniteElementSpace *f)
Associate a new FiniteElementSpace with the ParGridFunction.
Definition: pgridfunc.cpp:96
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
virtual void Mult(const Vector &x, Vector &y) const
Matrix multiplication: .
void AssembleSharedFaces(int skip_zeros=1)
virtual const SparseMatrix * GetRestrictionMatrix() const
Get the R matrix which restricts a local dof vector to true dof vector.
Definition: pfespace.hpp:379
void DeleteAll()
Delete the whole array.
Definition: array.hpp:854
virtual void AssembleDiagonal(Vector &diag) const
Computes the diagonal entries into diag. Typically, this operation only makes sense for linear Operat...
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
virtual void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)=0
virtual void Update(FiniteElementSpace *nfes=NULL)
Update the FiniteElementSpace and delete all data associated with the old one.
double * GetData()
Return the element data, i.e. the array A.
Definition: sparsemat.hpp:232
void ReduceRHS(const Vector &b, Vector &b_r) const
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P_test^t A P_trial.
DofTransformation * GetFaceNbrElementVDofs(int i, Array< int > &vdofs) const
Definition: pfespace.cpp:1457
void SetSize(int m, int n)
Definition: array.hpp:375
ParFiniteElementSpace * pfes
Points to the same object as fes.
void HypreStealOwnership(HypreParMatrix &A_hyp, SparseMatrix &A_diag)
Make A_hyp steal ownership of its diagonal part A_diag.
Definition: hypre.cpp:2744
void EliminateBC(const OperatorHandle &A_e, const Array< int > &ess_dof_list, const Vector &X, Vector &B) const
Eliminate essential dofs from the solution X into the r.h.s. B.
Definition: handle.cpp:340
virtual void FormRectangularSystemOperator(const Array< int > &trial_tdof_list, const Array< int > &test_tdof_list, OperatorHandle &A)=0
Data type sparse matrix.
Definition: sparsemat.hpp:50
virtual void AssembleDiagonal(Vector &diag) const
Assemble the diagonal of the bilinear form into diag. Note that diag is a true-dof Vector...
StaticCondensation * static_cond
Owned.
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:759
void Clear()
Definition: table.cpp:380
SparseMatrix * mat
Sparse matrix to be associated with the form. Owned.
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:154
double b
Definition: lissajous.cpp:42
const HYPRE_BigInt * GetFaceNbrGlobalDofMap()
Definition: pfespace.hpp:389
void EliminateRowsCols(OperatorHandle &A, const Array< int > &ess_dof_list)
Reset the OperatorHandle to be the eliminated part of A after elimination of the essential dofs ess_d...
Definition: handle.cpp:252
FaceElementTransformations * GetSharedFaceTransformations(int sf, bool fill2=true)
Definition: pmesh.cpp:3080
ParMesh * GetParMesh() const
Definition: pfespace.hpp:273
void EliminateVDofsInRHS(const Array< int > &vdofs, const Vector &x, Vector &b)
HypreParMatrix & GetParallelMatrix()
Return the parallel Schur complement matrix.
Definition: staticcond.hpp:159
void AbsMultTranspose(double a, const Vector &x, double b, Vector &y) const
Computes y = a * |At| * x + b * y, using entry-wise absolute values of the transpose of the matrix A...
Definition: hypre.cpp:1898
void Assemble(int skip_zeros=1)
Assemble the local matrix.
void TrueAddMult(const Vector &x, Vector &y, const double a=1.0) const
Compute y += a (P^t A P) x, where x and y are vectors on the true dofs.
Table * GetFaceToElementTable() const
Definition: mesh.cpp:6480
virtual void FormRectangularSystemMatrix(const Array< int > &trial_tdof_list, const Array< int > &test_tdof_list, OperatorHandle &A)
Return in A a parallel (on truedofs) version of this operator.
virtual void Update(FiniteElementSpace *nfes=NULL)
Update the FiniteElementSpace and delete all data associated with the old one.
ParFiniteElementSpace * test_pfes
Points to the same object as test_fes.
SparseMatrix * mat
Owned.
Set the diagonal value to one.
Definition: operator.hpp:50
void RAP(const DenseMatrix &A, const DenseMatrix &P, DenseMatrix &RAP)
Definition: densemat.cpp:3232
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:281
bool Finalized() const
Returns whether or not CSR format has been finalized.
Definition: sparsemat.hpp:552
void AddSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Definition: sparsemat.cpp:2722
Dynamic 2D array using row-major layout.
Definition: array.hpp:354
void SetSubVectorComplement(const Array< int > &dofs, const double val)
Set all vector entries NOT in the dofs Array to the given val.
Definition: vector.cpp:740
void MakeRectangularBlockDiag(MPI_Comm comm, HYPRE_BigInt glob_num_rows, HYPRE_BigInt glob_num_cols, HYPRE_BigInt *row_starts, HYPRE_BigInt *col_starts, SparseMatrix *diag)
Reset the OperatorHandle to hold a parallel rectangular block-diagonal matrix using the currently set...
Definition: handle.cpp:91
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:702
void ReduceSystem(Vector &x, Vector &b, Vector &X, Vector &B, int copy_interior=0) const
Set the reduced solution X and r.h.s B vectors from the full linear system solution x and r...
Definition: staticcond.cpp:416
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for the i&#39;th element. The returned indices are offsets into an ...
Definition: fespace.cpp:281
SparseMatrix * mat_e
Sparse Matrix used to store the eliminations from the b.c. Owned. .
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:161
ParFiniteElementSpace * domain_fes
Points to the same object as trial_fes.
bool IsIdentityProlongation(const Operator *P)
Definition: operator.hpp:720
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:712
virtual void MultTranspose(const Vector &x, Vector &y) const
Multiply a vector with the transposed matrix. y = At * x.
Definition: sparsemat.cpp:908
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:555
int GetNSharedFaces() const
Return the number of shared faces (3D), edges (2D), vertices (1D)
Definition: pmesh.cpp:3196
void MakeSquareBlockDiag(MPI_Comm comm, HYPRE_BigInt glob_size, HYPRE_BigInt *row_starts, SparseMatrix *diag)
Reset the OperatorHandle to hold a parallel square block-diagonal matrix using the currently set type...
Definition: handle.cpp:60
virtual const FiniteElement * GetFE(int i) const
Definition: pfespace.cpp:534
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:219
HYPRE_BigInt GlobalVSize() const
Definition: pfespace.hpp:279
MixedBilinearFormExtension * ext
ParGridFunction Xaux
Auxiliary objects used in TrueAddMult().
ParFiniteElementSpace * ParFESpace() const
Definition: pgridfunc.hpp:113
void Finalize()
Finalize the construction of the Schur complement matrix.
Definition: staticcond.cpp:233
void ParallelEliminateEssentialBC(const Array< int > &bdr_attr_is_ess, HypreParMatrix &A, const HypreParVector &X, HypreParVector &B) const
Eliminate essential boundary DOFs from a parallel assembled system.
HYPRE_Int HYPRE_BigInt
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
Definition: table.cpp:413
Table * GetFaceToAllElementTable() const
Definition: pmesh.cpp:2974
void Clear()
Clear the OperatorHandle, deleting the held Operator (if owned), while leaving the type id unchanged...
Definition: handle.hpp:124
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
void Distribute(const Vector *tv)
Definition: pgridfunc.cpp:141
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
Reconstruct a solution vector x (e.g. a GridFunction) from the solution X of a constrained linear sys...
Definition: operator.cpp:148
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
HYPRE_BigInt * GetDofOffsets() const
Definition: pfespace.hpp:277
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1086
double a
Definition: lissajous.cpp:41
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Form the linear system matrix A, see FormLinearSystem() for details.
virtual const Operator * GetProlongationMatrix() const
The returned Operator is owned by the FiniteElementSpace.
Definition: fespace.hpp:593
Operator::Type Type() const
Get the currently set operator type id.
Definition: handle.hpp:99
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:248
void SetOperatorOwner(bool own=true)
Set the ownership flag for the held Operator.
Definition: handle.hpp:120
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get()...
Definition: handle.hpp:104
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
SparseMatrix * mat_e
Owned.
int Size() const
Returns the number of TYPE I elements.
Definition: table.hpp:92
Ordering::Type GetOrdering() const
Return the ordering method.
Definition: fespace.hpp:721
ParFiniteElementSpace * range_fes
Points to the same object as test_fes.
virtual void AssembleDiagonal(Vector &diag) const
Computes the diagonal entries into diag. Typically, this operation only makes sense for linear Operat...
Definition: operator.hpp:131
FiniteElementSpace * fes
FE space on which the form lives. Not owned.
Vector Xaux
Auxiliary vectors used in TrueAddMult(): L-, L-, and T-vector, resp.
void Copy(Array &copy) const
Create a copy of the internal array to the provided copy.
Definition: array.hpp:864
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:709
virtual void FormRectangularLinearSystem(const Array< int > &trial_tdof_list, const Array< int > &test_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B)
Form the parallel linear system A X = B, corresponding to this mixed bilinear form and the linear for...
Operator * Ptr() const
Access the underlying Operator pointer.
Definition: handle.hpp:87
ParFiniteElementSpace * trial_pfes
Points to the same object as trial_fes.
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
MPI_Comm GetComm() const
Definition: pfespace.hpp:269
OperatorHandle p_mat
Vector data type.
Definition: vector.hpp:58
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)=0
ID for class HypreParMatrix.
Definition: operator.hpp:287
Hybridization * hybridization
Owned.
int * GetI()
Definition: table.hpp:113
HypreParMatrix * Dof_TrueDof_Matrix() const
The true dof-to-dof interpolation matrix.
Definition: pfespace.hpp:317
void SetEssentialTrueDofs(const Array< int > &ess_tdof_list)
Determine and save internally essential reduced true dofs.
Definition: staticcond.hpp:128
Array< int > vdofs
HYPRE_Int MultTranspose(HypreParVector &x, HypreParVector &y, double alpha=1.0, double beta=0.0) const
Computes y = alpha * A^t * x + beta * y.
Definition: hypre.cpp:1872
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
void ComputeSolution(const Vector &b, const Vector &sol_r, Vector &sol) const
Class for parallel meshes.
Definition: pmesh.hpp:32
HYPRE_BigInt * GetTrueDofOffsets() const
Definition: pfespace.hpp:278
const Table & GetElementToDofTable() const
Return a reference to the internal Table that stores the lists of scalar dofs, for each mesh element...
Definition: fespace.hpp:1097
OperatorHandle p_mat
Matrix and eliminated matrix.
bool Conforming() const
Definition: fespace.hpp:561
void MakePtAP(OperatorHandle &A, OperatorHandle &P)
Reset the OperatorHandle to hold the product P^t A P.
Definition: handle.cpp:123
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28
void MakeRAP(OperatorHandle &Rt, OperatorHandle &A, OperatorHandle &P)
Reset the OperatorHandle to hold the product R A P, where R = Rt^t.
Definition: handle.cpp:161
void Reset(OpType *A, bool own_A=true)
Reset the OperatorHandle to the given OpType pointer, A.
Definition: handle.hpp:145
virtual void FormRectangularLinearSystem(const Array< int > &trial_tdof_list, const Array< int > &test_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B)=0