MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
pbilinearform.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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  if (Xaux.ParFESpace() != pfes)
351  {
352  Xaux.SetSpace(pfes);
353  Yaux.SetSpace(pfes);
355  }
356 
357  Xaux.Distribute(&x);
358  if (ext)
359  {
360  ext->Mult(Xaux, Yaux);
361  }
362  else
363  {
364  MFEM_VERIFY(interior_face_integs.Size() == 0,
365  "the case of interior face integrators is not"
366  " implemented");
367  mat->Mult(Xaux, Yaux);
368  }
370  y.Add(a,Ytmp);
371 }
372 
374  const Array<int> &ess_tdof_list, Vector &x, Vector &b,
375  OperatorHandle &A, Vector &X, Vector &B, int copy_interior)
376 {
377  if (ext)
378  {
379  ext->FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
380  return;
381  }
382 
383  // Finish the matrix assembly and perform BC elimination, storing the
384  // eliminated part of the matrix.
385  FormSystemMatrix(ess_tdof_list, A);
386 
387  const Operator &P = *pfes->GetProlongationMatrix();
388  const SparseMatrix &R = *pfes->GetRestrictionMatrix();
389 
390  // Transform the system and perform the elimination in B, based on the
391  // essential BC values from x. Restrict the BC part of x in X, and set the
392  // non-BC part to zero. Since there is no good initial guess for the Lagrange
393  // multipliers, set X = 0.0 for hybridization.
394  if (static_cond)
395  {
396  // Schur complement reduction to the exposed dofs
397  static_cond->ReduceSystem(x, b, X, B, copy_interior);
398  }
399  else if (hybridization)
400  {
401  // Reduction to the Lagrange multipliers system
402  HypreParVector true_X(pfes), true_B(pfes);
403  P.MultTranspose(b, true_B);
404  R.Mult(x, true_X);
405  p_mat.EliminateBC(p_mat_e, ess_tdof_list, true_X, true_B);
406  R.MultTranspose(true_B, b);
407  hybridization->ReduceRHS(true_B, B);
408  X.SetSize(B.Size());
409  X = 0.0;
410  }
411  else
412  {
413  // Variational restriction with P
414  X.SetSize(P.Width());
415  B.SetSize(X.Size());
416  P.MultTranspose(b, B);
417  R.Mult(x, X);
418  p_mat.EliminateBC(p_mat_e, ess_tdof_list, X, B);
419  if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
420  }
421 }
422 
424  const Array<int> &vdofs, const Vector &x, Vector &b)
425 {
426  p_mat.EliminateBC(p_mat_e, vdofs, x, b);
427 }
428 
430  OperatorHandle &A)
431 {
432  if (ext)
433  {
434  ext->FormSystemMatrix(ess_tdof_list, A);
435  return;
436  }
437 
438  // Finish the matrix assembly and perform BC elimination, storing the
439  // eliminated part of the matrix.
440  if (static_cond)
441  {
443  {
444  static_cond->SetEssentialTrueDofs(ess_tdof_list);
447  }
449  }
450  else
451  {
452  if (mat)
453  {
454  const int remove_zeros = 0;
455  Finalize(remove_zeros);
456  MFEM_VERIFY(p_mat.Ptr() == NULL && p_mat_e.Ptr() == NULL,
457  "The ParBilinearForm must be updated with Update() before "
458  "re-assembling the ParBilinearForm.");
460  delete mat;
461  mat = NULL;
462  delete mat_e;
463  mat_e = NULL;
464  p_mat_e.EliminateRowsCols(p_mat, ess_tdof_list);
465  }
466  if (hybridization)
467  {
469  }
470  else
471  {
472  A = p_mat;
473  }
474  }
475 }
476 
478  const Vector &X, const Vector &b, Vector &x)
479 {
480  if (ext)
481  {
482  ext->RecoverFEMSolution(X, b, x);
483  return;
484  }
485 
486  const Operator &P = *pfes->GetProlongationMatrix();
487 
488  if (static_cond)
489  {
490  // Private dofs back solve
491  static_cond->ComputeSolution(b, X, x);
492  }
493  else if (hybridization)
494  {
495  // Primal unknowns recovery
496  HypreParVector true_X(pfes), true_B(pfes);
497  P.MultTranspose(b, true_B);
498  const SparseMatrix &R = *pfes->GetRestrictionMatrix();
499  R.Mult(x, true_X); // get essential b.c. from x
500  hybridization->ComputeSolution(true_B, X, true_X);
501  x.SetSize(P.Height());
502  P.Mult(true_X, x);
503  }
504  else
505  {
506  // Apply conforming prolongation
508  P.Mult(X, x);
509  }
510 }
511 
513 {
514  BilinearForm::Update(nfes);
515 
516  if (nfes)
517  {
518  pfes = dynamic_cast<ParFiniteElementSpace *>(nfes);
519  MFEM_VERIFY(pfes != NULL, "nfes must be a ParFiniteElementSpace!");
520  }
521 
522  p_mat.Clear();
523  p_mat_e.Clear();
524 }
525 
526 
528 {
529  // construct the block-diagonal matrix A
530  HypreParMatrix *A =
536  mat);
537 
540 
541  delete A;
542 
543  return rap;
544 }
545 
547 {
548  // construct the rectangular block-diagonal matrix dA
549  OperatorHandle dA(A.Type());
555  mat);
556 
557  OperatorHandle P_test(A.Type()), P_trial(A.Type());
558 
559  // TODO - construct the Dof_TrueDof_Matrix directly in the required format.
561  P_trial.ConvertFrom(trial_pfes->Dof_TrueDof_Matrix());
562 
563  A.MakeRAP(P_test, dA, P_trial);
564 }
565 
566 /// Compute y += a (P^t A P) x, where x and y are vectors on the true dofs
568  const double a) const
569 {
570  if (Xaux.ParFESpace() != trial_pfes)
571  {
574  }
575 
576  Xaux.Distribute(&x);
577  mat->Mult(Xaux, Yaux);
579 }
580 
582  const Array<int>
583  &trial_tdof_list,
584  const Array<int> &test_tdof_list,
585  OperatorHandle &A)
586 {
587  if (ext)
588  {
589  ext->FormRectangularSystemOperator(trial_tdof_list, test_tdof_list, A);
590  return;
591  }
592 
593  if (mat)
594  {
595  Finalize();
597  delete mat;
598  mat = NULL;
599  delete mat_e;
600  mat_e = NULL;
601  HypreParMatrix *temp =
602  p_mat.As<HypreParMatrix>()->EliminateCols(trial_tdof_list);
603  p_mat.As<HypreParMatrix>()->EliminateRows(test_tdof_list);
604  p_mat_e.Reset(temp, true);
605  }
606 
607  A = p_mat;
608 }
609 
611  const Array<int>
612  &trial_tdof_list,
613  const Array<int> &test_tdof_list, Vector &x,
614  Vector &b, OperatorHandle &A, Vector &X,
615  Vector &B)
616 {
617  if (ext)
618  {
619  ext->FormRectangularLinearSystem(trial_tdof_list, test_tdof_list,
620  x, b, A, X, B);
621  return;
622  }
623 
624  FormRectangularSystemMatrix(trial_tdof_list, test_tdof_list, A);
625 
626  const Operator *test_P = test_pfes->GetProlongationMatrix();
627  const SparseMatrix *trial_R = trial_pfes->GetRestrictionMatrix();
628 
631  test_P->MultTranspose(b, B);
632  trial_R->Mult(x, X);
633 
634  p_mat_e.As<HypreParMatrix>()->Mult(-1.0, X, 1.0, B);
635  B.SetSubVector(test_tdof_list, 0.0);
636 }
637 
639 {
640  MFEM_ASSERT(mat, "Matrix is not assembled");
641  MFEM_ASSERT(mat->Finalized(), "Matrix is not finalized");
644  HypreParMatrix* RAP = P->LeftDiagMult(*RA, range_fes->GetTrueDofOffsets());
645  delete RA;
646  return RAP;
647 }
648 
650 {
651  // construct the rectangular block-diagonal matrix dA
652  OperatorHandle dA(A.Type());
658  mat);
659 
661  OperatorHandle R_test_transpose(A.Type());
662  R_test_transpose.MakeRectangularBlockDiag(range_fes->GetComm(),
667  Rt);
668 
669  // TODO - construct the Dof_TrueDof_Matrix directly in the required format.
670  OperatorHandle P_trial(A.Type());
672 
673  A.MakeRAP(R_test_transpose, dA, P_trial);
674  delete Rt;
675 }
676 
678 {
679  if (ext)
680  {
681  Array<int> empty;
682  ext->FormRectangularSystemOperator(empty, empty, A);
683  return;
684  }
685 
686  mfem_error("not implemented!");
687 }
688 
690 const
691 {
692  MFEM_VERIFY(mat->Finalized(), "Local matrix needs to be finalized for "
693  "GetParBlocks");
694 
696 
697  blocks.SetSize(range_fes->GetVDim(), domain_fes->GetVDim());
698 
699  RLP->GetBlocks(blocks,
702 
703  delete RLP;
704 }
705 
706 }
707 
708 #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.
Ordering::Type GetOrdering() const
Return the ordering method.
Definition: fespace.hpp:599
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
void EliminateRowsCols(const Array< int > &rows_cols, const HypreParVector &X, HypreParVector &B)
Definition: hypre.cpp:2272
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:573
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:587
virtual int NumNonZeroElems() const
Returns the number of the nonzero elements in the matrix.
Definition: sparsemat.cpp:1584
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Definition: pfespace.cpp:1032
void ConvertFrom(OperatorHandle &A)
Convert the given OperatorHandle A to the currently set type id.
Definition: handle.cpp:200
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get()...
Definition: handle.hpp:104
int * GetJ()
Definition: table.hpp:114
virtual void AssembleDiagonal(Vector &diag) const
Assemble the diagonal of the bilinear form into diag. Note that diag is a true-dof Vector...
ParMesh * GetParMesh() const
Definition: pfespace.hpp:277
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
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(...
virtual const Operator * GetProlongationMatrix() const
The returned Operator is owned by the FiniteElementSpace.
Definition: pfespace.cpp:1153
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:513
int TrueVSize() const
Obsolete, kept for backward compatibility.
Definition: pfespace.hpp:434
ParGridFunction Xaux
Auxiliary objects used in TrueAddMult().
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:73
HypreParMatrix & GetParallelMatrix()
Return the parallel hybridized matrix.
int GetNSharedFaces() const
Return the number of shared faces (3D), edges (2D), vertices (1D)
Definition: pmesh.cpp:3127
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
Operator::Type Type() const
Get the currently set operator type id.
Definition: handle.hpp:99
void ParallelRAP(SparseMatrix &loc_A, OperatorHandle &A, bool steal_loc_A=false)
Compute parallel RAP operator and store it in A as a HypreParMatrix.
void Copy(Array &copy) const
Create a copy of the internal array to the provided copy.
Definition: array.hpp:856
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:285
A specialized ElementTransformation class representing a face and its two neighboring elements...
Definition: eltrans.hpp:480
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
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
int Size() const
Returns the size of the vector.
Definition: vector.hpp:200
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 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
int * GetI()
Return the array I.
Definition: sparsemat.hpp:222
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:923
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 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:87
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
void AssembleSharedFaces(int skip_zeros=1)
void ComputeSolution(const Vector &b, const Vector &sc_sol, Vector &sol) const
Definition: staticcond.cpp:476
void DeleteAll()
Delete the whole array.
Definition: array.hpp:846
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
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:94
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
virtual int GetTrueVSize() const
Return the number of local vector true dofs.
Definition: pfespace.hpp:289
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P_test^t A P_trial.
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:1894
void SetSize(int m, int n)
Definition: array.hpp:370
HYPRE_BigInt GlobalVSize() const
Definition: pfespace.hpp:283
ParFiniteElementSpace * pfes
Points to the same object as fes.
const FiniteElement * GetFaceNbrFE(int i) const
Definition: pfespace.cpp:1517
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: fespace.cpp:281
void HypreStealOwnership(HypreParMatrix &A_hyp, SparseMatrix &A_diag)
Make A_hyp steal ownership of its diagonal part A_diag.
Definition: hypre.cpp:2740
MPI_Comm GetComm() const
Definition: pfespace.hpp:273
virtual const FiniteElement * GetFE(int i) const
Definition: pfespace.cpp:535
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
Computes the diagonal entries into diag. Typically, this operation only makes sense for linear Operat...
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:67
StaticCondensation * static_cond
Owned.
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:751
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:441
ParFiniteElementSpace * ParFESpace() const
Return the parallel FE space associated with the ParBilinearForm.
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:399
void ComputeSolution(const Vector &b, const Vector &sol_r, Vector &sol) const
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:3011
void EliminateVDofsInRHS(const Array< int > &vdofs, const Vector &x, Vector &b)
HypreParMatrix * Dof_TrueDof_Matrix() const
The true dof-to-dof interpolation matrix.
Definition: pfespace.hpp:321
HypreParMatrix & GetParallelMatrix()
Return the parallel Schur complement matrix.
Definition: staticcond.hpp:159
virtual const Operator * GetProlongationMatrix() const
The returned Operator is owned by the FiniteElementSpace.
Definition: fespace.hpp:479
void Assemble(int skip_zeros=1)
Assemble the local matrix.
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.
HYPRE_BigInt GetMyDofOffset() const
Definition: pfespace.cpp:1143
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:590
bool Conforming() const
Definition: fespace.hpp:447
ParFiniteElementSpace * test_pfes
Points to the same object as test_fes.
int Size() const
Returns the number of TYPE I elements.
Definition: table.hpp:92
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:581
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:3213
HYPRE_BigInt * GetDofOffsets() const
Definition: pfespace.hpp:281
void AddSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Definition: sparsemat.cpp:2718
Operator * Ptr() const
Access the underlying Operator pointer.
Definition: handle.hpp:87
Dynamic 2D array using row-major layout.
Definition: array.hpp:351
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:708
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
bool Finalized() const
Returns whether or not CSR format has been finalized.
Definition: sparsemat.hpp:551
HypreParMatrix * ParallelAssemble() const
Returns the matrix &quot;assembled&quot; on the true dofs.
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.
HYPRE_BigInt * GetTrueDofOffsets() const
Definition: pfespace.hpp:282
bool IsIdentityProlongation(const Operator *P)
Definition: operator.hpp:693
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
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
MixedBilinearFormExtension * ext
ParGridFunction Xaux
Auxiliary objects used in TrueAddMult().
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Definition: sparsemat.cpp:729
void Finalize()
Finalize the construction of the Schur complement matrix.
Definition: staticcond.cpp:233
HYPRE_Int HYPRE_BigInt
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
Definition: table.cpp:413
void Clear()
Clear the OperatorHandle, deleting the held Operator (if owned), while leaving the type id unchanged...
Definition: handle.hpp:124
void Distribute(const Vector *tv)
Definition: pgridfunc.cpp:132
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:1868
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:85
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
Table * GetFaceToElementTable() const
Definition: mesh.cpp:6114
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.
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:252
void SetOperatorOwner(bool own=true)
Set the ownership flag for the held Operator.
Definition: handle.hpp:120
virtual void AssembleDiagonal(Vector &diag) const
Computes the diagonal entries into diag. Typically, this operation only makes sense for linear Operat...
Definition: operator.hpp:108
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
virtual const SparseMatrix * GetRestrictionMatrix() const
Get the R matrix which restricts a local dof vector to true dof vector.
Definition: pfespace.hpp:389
SparseMatrix * mat_e
Owned.
ParFiniteElementSpace * range_fes
Points to the same object as test_fes.
bool HasEliminatedBC() const
Return true if essential boundary conditions have been eliminated from the Schur complement matrix...
Definition: staticcond.hpp:142
void MultTranspose(const Vector &x, Vector &y) const
Multiply a vector with the transposed matrix. y = At * x.
Definition: sparsemat.cpp:904
FiniteElementSpace * fes
FE space on which the form lives. Not owned.
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.
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...
ParFiniteElementSpace * trial_pfes
Points to the same object as trial_fes.
int GetFaceNbrVSize() const
Definition: pfespace.hpp:394
OperatorHandle p_mat
Vector data type.
Definition: vector.hpp:60
void GetParBlocks(Array2D< HypreParMatrix * > &blocks) const
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)=0
ID for class HypreParMatrix.
Definition: operator.hpp:260
Hybridization * hybridization
Owned.
void ReduceRHS(const Vector &b, Vector &b_r) const
int * GetI()
Definition: table.hpp:113
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:750
void SetEssentialTrueDofs(const Array< int > &ess_tdof_list)
Determine and save internally essential reduced true dofs.
Definition: staticcond.hpp:128
Table * GetFaceToAllElementTable() const
Definition: pmesh.cpp:2905
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
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.
Class for parallel meshes.
Definition: pmesh.hpp:32
OperatorHandle p_mat
Matrix and eliminated matrix.
virtual void Mult(const Vector &x, Vector &y) const
Matrix multiplication: .
void MakePtAP(OperatorHandle &A, OperatorHandle &P)
Reset the OperatorHandle to hold the product P^t A P.
Definition: handle.cpp:123
ParGridFunction Yaux
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
DofTransformation * GetFaceNbrElementVDofs(int i, Array< int > &vdofs) const
Definition: pfespace.cpp:1464
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
ParFiniteElementSpace * ParFESpace() const
Definition: pgridfunc.hpp:113