MFEM  v3.2
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, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
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 = new 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 {
126  if (m == NULL) { return NULL; }
127 
128  MFEM_VERIFY(m->Finalized(), "local matrix needs to be finalized for "
129  "ParallelAssemble");
130 
131  HypreParMatrix *A;
132  if (fbfi.Size() == 0)
133  {
134  // construct a parallel block-diagonal wrapper matrix A based on m
135  A = new HypreParMatrix(pfes->GetComm(),
136  pfes->GlobalVSize(), pfes->GetDofOffsets(), m);
137  }
138  else
139  {
140  // handle the case when 'm' contains offdiagonal
141  int lvsize = pfes->GetVSize();
142  const HYPRE_Int *face_nbr_glob_ldof = pfes->GetFaceNbrGlobalDofMap();
143  HYPRE_Int ldof_offset = pfes->GetMyDofOffset();
144 
145  Array<HYPRE_Int> glob_J(m->NumNonZeroElems());
146  int *J = m->GetJ();
147  for (int i = 0; i < glob_J.Size(); i++)
148  {
149  if (J[i] < lvsize)
150  {
151  glob_J[i] = J[i] + ldof_offset;
152  }
153  else
154  {
155  glob_J[i] = face_nbr_glob_ldof[J[i] - lvsize];
156  }
157  }
158 
159  A = new HypreParMatrix(pfes->GetComm(), lvsize, pfes->GlobalVSize(),
160  pfes->GlobalVSize(), m->GetI(), glob_J,
161  m->GetData(), pfes->GetDofOffsets(),
162  pfes->GetDofOffsets());
163  }
164 
166 
167  delete A;
168 
169  return rap;
170 }
171 
173 {
174  ParMesh *pmesh = pfes->GetParMesh();
176  Array<int> vdofs1, vdofs2, vdofs_all;
178 
179  int nfaces = pmesh->GetNSharedFaces();
180  for (int i = 0; i < nfaces; i++)
181  {
182  T = pmesh->GetSharedFaceTransformations(i);
183  pfes->GetElementVDofs(T->Elem1No, vdofs1);
184  pfes->GetFaceNbrElementVDofs(T->Elem2No, vdofs2);
185  vdofs1.Copy(vdofs_all);
186  for (int j = 0; j < vdofs2.Size(); j++)
187  {
188  vdofs2[j] += height;
189  }
190  vdofs_all.Append(vdofs2);
191  for (int k = 0; k < fbfi.Size(); k++)
192  {
193  fbfi[k]->AssembleFaceMatrix(*pfes->GetFE(T->Elem1No),
194  *pfes->GetFaceNbrFE(T->Elem2No),
195  *T, elemmat);
196  if (keep_nbr_block)
197  {
198  mat->AddSubMatrix(vdofs_all, vdofs_all, elemmat, skip_zeros);
199  }
200  else
201  {
202  mat->AddSubMatrix(vdofs1, vdofs_all, elemmat, skip_zeros);
203  }
204  }
205  }
206 }
207 
208 void ParBilinearForm::Assemble(int skip_zeros)
209 {
210  if (mat == NULL && fbfi.Size() > 0)
211  {
213  pAllocMat();
214  }
215 
216  BilinearForm::Assemble(skip_zeros);
217 
218  if (fbfi.Size() > 0)
219  {
220  AssembleSharedFaces(skip_zeros);
221  }
222 }
223 
224 void ParBilinearForm
226  HypreParMatrix &A, const HypreParVector &X,
227  HypreParVector &B) const
228 {
229  Array<int> dof_list;
230 
231  pfes->GetEssentialTrueDofs(bdr_attr_is_ess, dof_list);
232 
233  // do the parallel elimination
234  A.EliminateRowsCols(dof_list, X, B);
235 }
236 
239  HypreParMatrix &A) const
240 {
241  Array<int> dof_list;
242 
243  pfes->GetEssentialTrueDofs(bdr_attr_is_ess, dof_list);
244 
245  return A.EliminateRowsCols(dof_list);
246 }
247 
248 void ParBilinearForm::TrueAddMult(const Vector &x, Vector &y, const double a)
249 const
250 {
251  MFEM_VERIFY(fbfi.Size() == 0, "the case of interior face integrators is not"
252  " implemented");
253 
254  if (X.ParFESpace() != pfes)
255  {
256  X.SetSpace(pfes);
257  Y.SetSpace(pfes);
258  }
259 
260  X.Distribute(&x);
261  mat->Mult(X, Y);
262  pfes->Dof_TrueDof_Matrix()->MultTranspose(a, Y, 1.0, y);
263 }
264 
266  Array<int> &ess_tdof_list, Vector &x, Vector &b,
267  HypreParMatrix &A, Vector &X, Vector &B, int copy_interior)
268 {
270  const SparseMatrix &R = *pfes->GetRestrictionMatrix();
271  Array<int> ess_rtdof_list;
272 
273  // Finish the matrix assembly and perform BC elimination, storing the
274  // eliminated part of the matrix.
275  if (static_cond)
276  {
277  static_cond->ConvertListToReducedTrueDofs(ess_tdof_list, ess_rtdof_list);
279  {
281  static_cond->EliminateReducedTrueDofs(ess_rtdof_list, 0);
282  }
283  }
284  else if (mat)
285  {
286  MFEM_VERIFY(p_mat == NULL && p_mat_e == NULL,
287  "The ParBilinearForm must be updated with Update() before "
288  "re-assembling the ParBilinearForm.");
289  Finalize();
291  delete mat;
292  mat = NULL;
293  delete mat_e;
294  mat_e = NULL;
295  p_mat_e = p_mat->EliminateRowsCols(ess_tdof_list);
296  }
297 
298  // Transform the system and perform the elimination in B, based on the
299  // essential BC values from x. Restrict the BC part of x in X, and set the
300  // non-BC part to zero. Since there is no good initial guess for the Lagrange
301  // multipliers, set X = 0.0 for hybridization.
302  if (static_cond)
303  {
304  // Schur complement reduction to the exposed dofs
305  static_cond->ReduceRHS(b, B);
309  ess_rtdof_list, X, B);
310  if (!copy_interior) { X.SetSubVectorComplement(ess_rtdof_list, 0.0); }
312  }
313  else if (hybridization)
314  {
315  // Reduction to the Lagrange multipliers system
316  HypreParVector true_X(pfes), true_B(pfes);
317  P.MultTranspose(b, true_B);
318  R.Mult(x, true_X);
319  EliminateBC(*p_mat, *p_mat_e, ess_tdof_list, true_X, true_B);
320  R.MultTranspose(true_B, b);
321  hybridization->ReduceRHS(true_B, B);
322  X.SetSize(B.Size());
323  X = 0.0;
325  }
326  else
327  {
328  // Variational restriction with P
329  X.SetSize(pfes->TrueVSize());
330  B.SetSize(X.Size());
331  P.MultTranspose(b, B);
332  R.Mult(x, X);
333  EliminateBC(*p_mat, *p_mat_e, ess_tdof_list, X, B);
334  if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
335  A.MakeRef(*p_mat);
336  }
337 }
338 
340  const Vector &X, const Vector &b, Vector &x)
341 {
343 
344  if (static_cond)
345  {
346  // Private dofs back solve
347  static_cond->ComputeSolution(b, X, x);
348  }
349  else if (hybridization)
350  {
351  // Primal unknowns recovery
352  HypreParVector true_X(pfes), true_B(pfes);
353  P.MultTranspose(b, true_B);
354  const SparseMatrix &R = *pfes->GetRestrictionMatrix();
355  R.Mult(x, true_X); // get essential b.c. from x
356  hybridization->ComputeSolution(true_B, X, true_X);
357  x.SetSize(P.Height());
358  P.Mult(true_X, x);
359  }
360  else
361  {
362  // Apply conforming prolongation
363  x.SetSize(P.Height());
364  P.Mult(X, x);
365  }
366 }
367 
369 {
370  BilinearForm::Update(nfes);
371 
372  if (nfes)
373  {
374  pfes = dynamic_cast<ParFiniteElementSpace *>(nfes);
375  MFEM_VERIFY(pfes != NULL, "nfes must be a ParFiniteElementSpace!");
376  }
377 
378  delete p_mat;
379  delete p_mat_e;
380  p_mat = p_mat_e = NULL;
381 }
382 
383 
385 {
386  MFEM_ASSERT(mat, "matrix is not assembled");
387  MFEM_ASSERT(mat->Finalized(), "matrix is not finalized");
391  delete RA;
392  return RAP;
393 }
394 
396 const
397 {
398  MFEM_VERIFY(mat->Finalized(), "local matrix needs to be finalized for "
399  "GetParBlocks");
400 
402 
403  blocks.SetSize(range_fes->GetVDim(), domain_fes->GetVDim());
404 
405  RLP->GetBlocks(blocks,
408 
409  delete RLP;
410 }
411 
413 {
414  // construct the block-diagonal matrix A
415  HypreParMatrix *A =
421  mat);
422 
425 
426  delete A;
427 
428  return rap;
429 }
430 
433  const double a) const
434 {
435  if (X.ParFESpace() != trial_pfes)
436  {
439  }
440 
441  X.Distribute(&x);
442  mat->Mult(X, Y);
444 }
445 
446 }
447 
448 #endif
HypreParMatrix * p_mat_e
Ordering::Type GetOrdering() const
Return the ordering method.
Definition: fespace.hpp:173
int Size() const
Logical size of the array.
Definition: array.hpp:109
void EliminateRowsCols(const Array< int > &rows_cols, const HypreParVector &X, HypreParVector &B)
Definition: hypre.cpp:1260
int GetVSize() const
Definition: fespace.hpp:161
virtual int NumNonZeroElems() const
Returns the number of the nonzero elements in the matrix.
Definition: sparsemat.cpp:905
int * GetJ()
Definition: table.hpp:108
void SetSpace(ParFiniteElementSpace *f)
Definition: pgridfunc.cpp:71
Array< BilinearFormIntegrator * > fbfi
Set of interior face Integrators to be applied.
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
void ReduceRHS(const Vector &b, Vector &sc_b) const
Definition: staticcond.cpp:309
void MakeRef(const HypreParMatrix &master)
Make this HypreParMatrix a reference to &#39;master&#39;.
Definition: hypre.cpp:763
HypreParMatrix * LeftDiagMult(const SparseMatrix &D, HYPRE_Int *row_starts=NULL) const
Definition: hypre.cpp:1016
void SetSize(int s)
Resize the vector if the new size is different.
Definition: vector.hpp:263
HYPRE_Int MultTranspose(HypreParVector &x, HypreParVector &y, double alpha=1.0, double beta=0.0)
Computes y = alpha * A^t * x + beta * y.
Definition: hypre.cpp:1010
void GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: fespace.cpp:132
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:451
HypreParMatrix & GetParallelMatrix()
Return the parallel hybridized matrix.
int GetNSharedFaces() const
Return the number of shared faces (3D), edges (2D), vertices (1D)
Definition: pmesh.cpp:1547
void SetDims(int rows, int nnz)
Definition: table.cpp:132
void Copy(Array &copy) const
Create a copy of the current array.
Definition: array.hpp:160
HYPRE_Int * GetDofOffsets()
Definition: pfespace.hpp:166
Data type dense matrix using column-major storage.
Definition: densemat.hpp:22
int Size() const
Returns the size of the vector.
Definition: vector.hpp:86
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 LoseData()
Call this if data has been stolen.
Definition: table.hpp:141
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual const SparseMatrix * GetRestrictionMatrix()
Get the R matrix which restricts a local dof vector to true dof vector.
Definition: pfespace.hpp:235
void AssembleSharedFaces(int skip_zeros=1)
void ComputeSolution(const Vector &b, const Vector &sc_sol, Vector &sol) const
Definition: staticcond.cpp:452
const HYPRE_Int * GetFaceNbrGlobalDofMap()
Definition: pfespace.hpp:245
HYPRE_Int Mult(HypreParVector &x, HypreParVector &y, double alpha=1.0, double beta=0.0)
Computes y = alpha * A * x + beta * y.
Definition: hypre.cpp:940
HYPRE_Int GetMyDofOffset() const
Definition: pfespace.cpp:571
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
virtual void Update(FiniteElementSpace *nfes=NULL)
HYPRE_Int * GetTrueDofOffsets()
Definition: pfespace.hpp:167
HypreParMatrix * p_mat
void EliminateBC(HypreParMatrix &A, HypreParMatrix &Ae, const Array< int > &ess_dof_list, const Vector &X, Vector &B)
Definition: hypre.cpp:1415
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
void GetFaceNbrElementVDofs(int i, Array< int > &vdofs) const
Definition: pfespace.cpp:801
void SetSize(int m, int n)
Definition: array.hpp:264
ParFiniteElementSpace * pfes
const FiniteElement * GetFaceNbrFE(int i) const
Definition: pfespace.cpp:835
Data type sparse matrix.
Definition: sparsemat.hpp:38
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows.
Definition: operator.hpp:35
StaticCondensation * static_cond
int Append(const T &el)
Append element to array, resize if necessary.
Definition: array.hpp:376
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:136
void Clear()
Definition: table.cpp:346
SparseMatrix * mat
Sparse matrix to be associated with the form.
void ComputeSolution(const Vector &b, const Vector &sol_r, Vector &sol) const
FaceElementTransformations * GetSharedFaceTransformations(int sf, bool fill2=true)
Definition: pmesh.cpp:1462
HypreParMatrix * RAP(HypreParMatrix *A, HypreParMatrix *P)
Returns the matrix P^t * A * P.
Definition: hypre.cpp:1365
HypreParMatrix & GetParallelMatrix()
Return the parallel Schur complement matrix.
Definition: staticcond.hpp:148
void Assemble(int skip_zeros=1)
Assemble the local matrix.
void Transpose(const Table &A, Table &At, int _ncols_A)
Transpose a Table.
Definition: table.cpp:391
virtual void Update(FiniteElementSpace *nfes=NULL)
ParFiniteElementSpace * test_pfes
int Size() const
Returns the number of TYPE I elements.
Definition: table.hpp:86
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:151
double * GetData() const
Return element data.
Definition: sparsemat.hpp:121
void AddSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Definition: sparsemat.cpp:1759
int * GetI() const
Return the array I.
Definition: sparsemat.hpp:117
void SetSubVectorComplement(const Array< int > &dofs, const double val)
Set all vector entries NOT in the &#39;dofs&#39; array to the given &#39;val&#39;.
Definition: vector.cpp:575
HypreParMatrix * Dof_TrueDof_Matrix()
The true dof-to-dof interpolation matrix.
Definition: pfespace.cpp:345
bool Finalized() const
Definition: sparsemat.hpp:261
HypreParMatrix * ParallelAssemble() const
Returns the matrix &quot;assembled&quot; on the true dofs.
SparseMatrix * mat_e
Matrix used to eliminate b.c.
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:58
ParFiniteElementSpace * domain_fes
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Definition: sparsemat.cpp:439
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list)
Definition: pfespace.cpp:479
void Finalize()
Finalize the construction of the Schur complement matrix.
Definition: staticcond.cpp:235
void Distribute(const Vector *tv)
Definition: pgridfunc.cpp:85
void EliminateReducedTrueDofs(const Array< int > &ess_rtdof_list, int keep_diagonal)
Eliminate the given reduced true dofs from the Schur complement matrix S.
Definition: staticcond.cpp:286
void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
HypreParMatrix & GetParallelMatrixElim()
Return the eliminated part of the parallel Schur complement matrix.
Definition: staticcond.hpp:151
Table * GetFaceToElementTable() const
Definition: mesh.cpp:3531
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
ParFiniteElementSpace * range_fes
bool HasEliminatedBC() const
Definition: staticcond.hpp:131
void FormLinearSystem(Array< int > &ess_tdof_list, Vector &x, Vector &b, HypreParMatrix &A, Vector &X, Vector &B, int copy_interior=0)
void MultTranspose(const Vector &x, Vector &y) const
Multiply a vector with the transposed matrix. y = At * x.
Definition: sparsemat.cpp:514
FiniteElementSpace * fes
FE space on which the form lives.
void ParallelEliminateEssentialBC(const Array< int > &bdr_attr_is_ess, HypreParMatrix &A, const HypreParVector &X, HypreParVector &B) const
const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement associated with i&#39;th element.
Definition: fespace.cpp:1113
DenseMatrix elemmat
ParFiniteElementSpace * trial_pfes
int GetFaceNbrVSize() const
Definition: pfespace.hpp:240
Vector data type.
Definition: vector.hpp:33
void GetParBlocks(Array2D< HypreParMatrix * > &blocks) const
Hybridization * hybridization
void ReduceRHS(const Vector &b, Vector &b_r) const
int * GetI()
Definition: table.hpp:107
const Table & GetElementToDofTable() const
Definition: fespace.hpp:275
void ReduceSolution(const Vector &sol, Vector &sc_sol) const
Definition: staticcond.cpp:388
Table * GetFaceToAllElementTable() const
Definition: pmesh.cpp:1374
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:150
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.
int * GetJ() const
Return the array J.
Definition: sparsemat.hpp:119
Class for parallel meshes.
Definition: pmesh.hpp:28
void ConvertListToReducedTrueDofs(const Array< int > &ess_tdof_list, Array< int > &ess_rtdof_list) const
Definition: staticcond.hpp:169
ParFiniteElementSpace * ParFESpace() const
Definition: pgridfunc.hpp:61