MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
hybridization.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 "hybridization.hpp"
13 #include "gridfunc.hpp"
14 
15 #ifdef MFEM_USE_MPI
16 #include "pfespace.hpp"
17 #endif
18 
19 #include <map>
20 
21 // uncomment next line for debugging: write C and P to file
22 // #define MFEM_DEBUG_HYBRIDIZATION_CP
23 #ifdef MFEM_DEBUG_HYBRIDIZATION_CP
24 #include <fstream>
25 #endif
26 
27 namespace mfem
28 {
29 
31  FiniteElementSpace *c_fespace)
32  : fes(fespace), c_fes(c_fespace), c_bfi(NULL), Ct(NULL), H(NULL),
33  Af_data(NULL), Af_ipiv(NULL)
34 {
35 #ifdef MFEM_USE_MPI
36  pC = P_pc = NULL;
38 #endif
39 }
40 
42 {
43 #ifdef MFEM_USE_MPI
44  delete P_pc;
45  delete pC;
46 #endif
47  delete [] Af_ipiv;
48  delete [] Af_data;
49  delete H;
50  delete Ct;
51  delete c_bfi;
52 }
53 
55 {
56  const int NE = fes->GetNE();
57  int num_hat_dofs = hat_offsets[NE];
58  Array<int> vdofs, c_vdofs;
59 
60  int c_num_face_nbr_dofs = 0;
61 #ifdef MFEM_USE_MPI
62  ParFiniteElementSpace *c_pfes = dynamic_cast<ParFiniteElementSpace*>(c_fes);
63  ParMesh *pmesh = c_pfes ? c_pfes->GetParMesh() : NULL;
64  HYPRE_BigInt num_shared_slave_faces = 0, glob_num_shared_slave_faces = 0;
65  if (c_pfes)
66  {
67  if (pmesh->Nonconforming())
68  {
69  const int dim = pmesh->Dimension();
70  const NCMesh::NCList &shared = pmesh->pncmesh->GetSharedList(dim-1);
71  num_shared_slave_faces = (HYPRE_BigInt) shared.slaves.Size();
72  MPI_Allreduce(&num_shared_slave_faces, &glob_num_shared_slave_faces, 1,
73  HYPRE_MPI_BIG_INT, MPI_SUM, pmesh->GetComm());
74  MFEM_ASSERT(glob_num_shared_slave_faces%2 == 0, "");
75  glob_num_shared_slave_faces /= 2;
76  if (glob_num_shared_slave_faces)
77  {
78  c_pfes->ExchangeFaceNbrData();
79  c_num_face_nbr_dofs = c_pfes->GetFaceNbrVSize();
80  }
81 #ifdef MFEM_DEBUG_HERE
82  MFEM_WARNING('[' << c_pfes->GetMyRank() <<
83  "] num_shared_slave_faces = " << num_shared_slave_faces
84  << ", glob_num_shared_slave_faces = "
85  << glob_num_shared_slave_faces
86  << "\n num_face_nbr_dofs = " << c_num_face_nbr_dofs
87  << ", num_shared_faces = " << pmesh->GetNSharedFaces());
88 #undef MFEM_DEBUG_HERE
89 #endif
90  }
91  }
92 #endif
93 
94  const int c_vsize = c_fes->GetVSize();
95  Ct = new SparseMatrix(num_hat_dofs, c_vsize + c_num_face_nbr_dofs);
96 
97  if (c_bfi)
98  {
99  const int skip_zeros = 1;
100  DenseMatrix elmat;
102  Mesh *mesh = fes->GetMesh();
103  int num_faces = mesh->GetNumFaces();
104  for (int i = 0; i < num_faces; i++)
105  {
106  FTr = mesh->GetInteriorFaceTransformations(i);
107  if (!FTr) { continue; }
108 
109  int o1 = hat_offsets[FTr->Elem1No];
110  int s1 = hat_offsets[FTr->Elem1No+1] - o1;
111  int o2 = hat_offsets[FTr->Elem2No];
112  int s2 = hat_offsets[FTr->Elem2No+1] - o2;
113  vdofs.SetSize(s1 + s2);
114  for (int j = 0; j < s1; j++)
115  {
116  vdofs[j] = o1 + j;
117  }
118  for (int j = 0; j < s2; j++)
119  {
120  vdofs[s1+j] = o2 + j;
121  }
122  c_fes->GetFaceVDofs(i, c_vdofs);
124  *fes->GetFE(FTr->Elem1No),
125  *fes->GetFE(FTr->Elem2No),
126  *FTr, elmat);
127  // zero-out small elements in elmat
128  elmat.Threshold(1e-12 * elmat.MaxMaxNorm());
129  Ct->AddSubMatrix(vdofs, c_vdofs, elmat, skip_zeros);
130  }
131 #ifdef MFEM_USE_MPI
132  if (pmesh)
133  {
134  // Assemble local contribution to Ct from shared faces
135  const int num_shared_faces = pmesh->GetNSharedFaces();
136  for (int i = 0; i < num_shared_faces; i++)
137  {
138  const int face_no = pmesh->GetSharedFace(i);
139  const bool ghost_sface = (face_no >= num_faces);
140  const FiniteElement *fe, *face_fe;
141  if (!ghost_sface)
142  {
143  FTr = pmesh->GetFaceElementTransformations(face_no);
144  MFEM_ASSERT(FTr->Elem2No < 0, "");
145  face_fe = c_fes->GetFaceElement(face_no);
146  c_fes->GetFaceVDofs(face_no, c_vdofs);
147  }
148  else
149  {
150  const int fill2 = false; // only need side "1" data
151  FTr = pmesh->GetSharedFaceTransformations(i, fill2);
152  face_fe = c_pfes->GetFaceNbrFaceFE(face_no);
153  c_pfes->GetFaceNbrFaceVDofs(face_no, c_vdofs);
155  for (int j = 0; j < c_vdofs.Size(); j++)
156  {
157  c_vdofs[j] += c_vsize;
158  }
159  }
160  int o1 = hat_offsets[FTr->Elem1No];
161  int s1 = hat_offsets[FTr->Elem1No+1] - o1;
162  vdofs.SetSize(s1);
163  for (int j = 0; j < s1; j++)
164  {
165  vdofs[j] = o1 + j;
166  }
167  fe = fes->GetFE(FTr->Elem1No);
168  c_bfi->AssembleFaceMatrix(*face_fe, *fe, *fe, *FTr, elmat);
169  // zero-out small elements in elmat
170  elmat.Threshold(1e-12 * elmat.MaxMaxNorm());
171  Ct->AddSubMatrix(vdofs, c_vdofs, elmat, skip_zeros);
172  }
173  if (glob_num_shared_slave_faces)
174  {
175  // Convert Ct to parallel and then transpose it:
176  Ct->Finalize(skip_zeros);
177  HYPRE_BigInt Ct_num_rows = Ct->Height();
178  Array<HYPRE_BigInt> Ct_rows;
179  Array<HYPRE_BigInt> *offsets[1] = { &Ct_rows };
180  pmesh->GenerateOffsets(1, &Ct_num_rows, offsets);
182  HYPRE_BigInt c_ldof_offset = c_pfes->GetMyDofOffset();
183  const HYPRE_BigInt *c_face_nbr_glob_ldof =
184  c_pfes->GetFaceNbrGlobalDofMap();
185  int *J = Ct->GetJ();
186  for (int i = 0; i < Ct_J.Size(); i++)
187  {
188  Ct_J[i] = J[i] < c_vsize ?
189  J[i] + c_ldof_offset :
190  c_face_nbr_glob_ldof[J[i] - c_vsize];
191  }
192  HypreParMatrix pCt(pmesh->GetComm(), Ct->Height(),
193  Ct_rows.Last(), c_pfes->GlobalVSize(),
194  Ct->GetI(), Ct_J.GetData(), Ct->GetData(),
195  Ct_rows, c_pfes->GetDofOffsets());
196  Ct_J.DeleteAll();
197  pC = pCt.Transpose();
198  }
199  if (pmesh->Nonconforming())
200  {
201  // TODO - Construct P_pc directly in the pH format
203  }
204  }
205 #endif
206  Ct->Finalize(skip_zeros);
207  }
208  else
209  {
210  // Check if c_fes is really needed here.
211  MFEM_ABORT("TODO: algebraic definition of C");
212  }
213 }
214 
215 void Hybridization::Init(const Array<int> &ess_tdof_list)
216 {
217  if (Ct) { return; }
218 
219  // count the number of dofs in the discontinuous version of fes:
220  const int NE = fes->GetNE();
221  Array<int> vdofs;
222  int num_hat_dofs = 0;
223  hat_offsets.SetSize(NE+1);
224  hat_offsets[0] = 0;
225  for (int i = 0; i < NE; i++)
226  {
227  fes->GetElementVDofs(i, vdofs);
228  num_hat_dofs += vdofs.Size();
229  hat_offsets[i+1] = num_hat_dofs;
230  }
231 
232  // Assemble the constraint matrix C
233  ConstructC();
234 
235 #ifdef MFEM_DEBUG_HYBRIDIZATION_CP
236  // Debug: write C and P to file
237  {
238  std::ofstream C_file("C_matrix.txt");
239  SparseMatrix *C = Transpose(*Ct);
240  C->PrintMatlab(C_file);
241  delete C;
242 
244  if (P)
245  {
246  std::ofstream P_file("P_matrix.txt");
247  P->PrintMatlab(P_file);
248  }
249  }
250 #endif
251 
252  // Define the "free" (0) and "essential" (1) hat_dofs.
253  // The "essential" hat_dofs are those that depend only on essential cdofs;
254  // all other hat_dofs are "free".
255  hat_dofs_marker.SetSize(num_hat_dofs);
256  Array<int> free_tdof_marker;
257 #ifdef MFEM_USE_MPI
258  ParFiniteElementSpace *pfes = dynamic_cast<ParFiniteElementSpace*>(fes);
259  free_tdof_marker.SetSize(pfes ? pfes->TrueVSize() :
261 #else
262  free_tdof_marker.SetSize(fes->GetConformingVSize());
263 #endif
264  free_tdof_marker = 1;
265  for (int i = 0; i < ess_tdof_list.Size(); i++)
266  {
267  free_tdof_marker[ess_tdof_list[i]] = 0;
268  }
269  Array<int> free_vdofs_marker;
270 #ifdef MFEM_USE_MPI
271  if (!pfes)
272  {
274  if (!cP)
275  {
276  free_vdofs_marker.MakeRef(free_tdof_marker);
277  }
278  else
279  {
280  free_vdofs_marker.SetSize(fes->GetVSize());
281  cP->BooleanMult(free_tdof_marker, free_vdofs_marker);
282  }
283  }
284  else
285  {
286  HypreParMatrix *P = pfes->Dof_TrueDof_Matrix();
287  free_vdofs_marker.SetSize(fes->GetVSize());
288  P->BooleanMult(1, free_tdof_marker, 0, free_vdofs_marker);
289  }
290 #else
292  if (!cP)
293  {
294  free_vdofs_marker.MakeRef(free_tdof_marker);
295  }
296  else
297  {
298  free_vdofs_marker.SetSize(fes->GetVSize());
299  cP->BooleanMult(free_tdof_marker, free_vdofs_marker);
300  }
301 #endif
302  for (int i = 0; i < NE; i++)
303  {
304  fes->GetElementVDofs(i, vdofs);
306  for (int j = 0; j < vdofs.Size(); j++)
307  {
308  hat_dofs_marker[hat_offsets[i]+j] = ! free_vdofs_marker[vdofs[j]];
309  }
310  }
311 #ifndef MFEM_DEBUG
312  // In DEBUG mode this array is used below.
313  free_tdof_marker.DeleteAll();
314 #endif
315  free_vdofs_marker.DeleteAll();
316  // Split the "free" (0) hat_dofs into "internal" (0) or "boundary" (-1).
317  // The "internal" hat_dofs are those "free" hat_dofs for which the
318  // corresponding column in C is zero; otherwise the free hat_dof is
319  // "boundary".
320  for (int i = 0; i < num_hat_dofs; i++)
321  {
322  // skip "essential" hat_dofs and empty rows in Ct
323  if (hat_dofs_marker[i] != 1 && Ct->RowSize(i) > 0)
324  {
325  hat_dofs_marker[i] = -1; // mark this hat_dof as "boundary"
326  }
327  }
328 
329  // Define Af_offsets and Af_f_offsets
330  Af_offsets.SetSize(NE+1);
331  Af_offsets[0] = 0;
332  Af_f_offsets.SetSize(NE+1);
333  Af_f_offsets[0] = 0;
334  // #define MFEM_DEBUG_HERE // uncomment to enable printing of hat dofs stats
335 #ifdef MFEM_DEBUG_HERE
336  int b_size = 0;
337 #endif
338  for (int i = 0; i < NE; i++)
339  {
340  int f_size = 0; // count the "free" hat_dofs in element i
341  for (int j = hat_offsets[i]; j < hat_offsets[i+1]; j++)
342  {
343  if (hat_dofs_marker[j] != 1) { f_size++; }
344 #ifdef MFEM_DEBUG_HERE
345  if (hat_dofs_marker[j] == -1) { b_size++; }
346 #endif
347  }
348  Af_offsets[i+1] = Af_offsets[i] + f_size*f_size;
349  Af_f_offsets[i+1] = Af_f_offsets[i] + f_size;
350  }
351 
352 #ifdef MFEM_DEBUG_HERE
353 #ifndef MFEM_USE_MPI
354  int myid = 0;
355 #else
356  int myid = pmesh ? pmesh->GetMyRank() : 0;
357 #endif
358  int i_size = Af_f_offsets[NE] - b_size;
359  int e_size = num_hat_dofs - (i_size + b_size);
360  mfem::out << "\nHybridization::Init:"
361  << " [" << myid << "] hat dofs - \"internal\": " << i_size
362  << ", \"boundary\": " << b_size
363  << ", \"essential\": " << e_size << '\n' << std::endl;
364 #undef MFEM_DEBUG_HERE
365 #endif
366 
367  Af_data = new double[Af_offsets[NE]];
368  Af_ipiv = new int[Af_f_offsets[NE]];
369 
370 #ifdef MFEM_DEBUG
371  // check that Ref = 0
372  const SparseMatrix *R = fes->GetRestrictionMatrix();
373  if (!R) { return; }
374  Array<int> vdof_marker(fes->GetVSize()); // 0 - f, 1 - e
375  vdof_marker = 0;
376  for (int i = 0; i < NE; i++)
377  {
378  fes->GetElementVDofs(i, vdofs);
380  for (int j = 0; j < vdofs.Size(); j++)
381  {
382  if (hat_dofs_marker[hat_offsets[i]+j] == 1) // "essential" hat dof
383  {
384  vdof_marker[vdofs[j]] = 1;
385  }
386  }
387  }
388  for (int tdof = 0; tdof < R->Height(); tdof++)
389  {
390  if (free_tdof_marker[tdof]) { continue; }
391 
392  const int ncols = R->RowSize(tdof);
393  const int *cols = R->GetRowColumns(tdof);
394  const double *vals = R->GetRowEntries(tdof);
395  for (int j = 0; j < ncols; j++)
396  {
397  if (std::abs(vals[j]) != 0.0 && vdof_marker[cols[j]] == 0)
398  {
399  MFEM_ABORT("Ref != 0");
400  }
401  }
402  }
403 #endif
404 }
405 
407  int el, Array<int> &i_dofs, Array<int> &b_dofs) const
408 {
409  // returns local indices in i_dofs and b_dofs
410  int h_start, h_end;
411 
412  h_start = hat_offsets[el];
413  h_end = hat_offsets[el+1];
414  i_dofs.Reserve(h_end-h_start);
415  i_dofs.SetSize(0);
416  b_dofs.Reserve(h_end-h_start);
417  b_dofs.SetSize(0);
418  for (int i = h_start; i < h_end; i++)
419  {
420  int mark = hat_dofs_marker[i];
421  if (mark == 0) { i_dofs.Append(i-h_start); }
422  else if (mark == -1) { b_dofs.Append(i-h_start); }
423  }
424 }
425 
426 void Hybridization::GetBDofs(int el, int &num_idofs, Array<int> &b_dofs) const
427 {
428  // returns global indices in b_dofs
429  const int h_start = hat_offsets[el];
430  const int h_end = hat_offsets[el+1];
431  b_dofs.Reserve(h_end-h_start);
432  b_dofs.SetSize(0);
433  num_idofs = 0;
434  for (int i = h_start; i < h_end; i++)
435  {
436  int mark = hat_dofs_marker[i];
437  if (mark == 0) { num_idofs++; }
438  else if (mark == -1) { b_dofs.Append(i); }
439  }
440 }
441 
443 {
444  Array<int> i_dofs, b_dofs;
445 
446  GetIBDofs(el, i_dofs, b_dofs);
447 
448  DenseMatrix A_ii(Af_data + Af_offsets[el], i_dofs.Size(), i_dofs.Size());
449  DenseMatrix A_ib(A_ii.Data() + i_dofs.Size()*i_dofs.Size(),
450  i_dofs.Size(), b_dofs.Size());
451  DenseMatrix A_bi(A_ib.Data() + i_dofs.Size()*b_dofs.Size(),
452  b_dofs.Size(), i_dofs.Size());
453  DenseMatrix A_bb(A_bi.Data() + b_dofs.Size()*i_dofs.Size(),
454  b_dofs.Size(), b_dofs.Size());
455 
456  for (int j = 0; j < i_dofs.Size(); j++)
457  {
458  int j_dof = i_dofs[j];
459  for (int i = 0; i < i_dofs.Size(); i++)
460  {
461  A_ii(i,j) = A(i_dofs[i],j_dof);
462  }
463  for (int i = 0; i < b_dofs.Size(); i++)
464  {
465  A_bi(i,j) = A(b_dofs[i],j_dof);
466  }
467  }
468  for (int j = 0; j < b_dofs.Size(); j++)
469  {
470  int j_dof = b_dofs[j];
471  for (int i = 0; i < i_dofs.Size(); i++)
472  {
473  A_ib(i,j) = A(i_dofs[i],j_dof);
474  }
475  for (int i = 0; i < b_dofs.Size(); i++)
476  {
477  A_bb(i,j) = A(b_dofs[i],j_dof);
478  }
479  }
480 }
481 
483 {
484  // Not tested.
485 #ifdef MFEM_DEBUG
486  Array<int> vdofs, bvdofs;
487  fes->GetBdrElementVDofs(bdr_el, bvdofs);
488 #endif
489 
490  int el;
491  DenseMatrix B(A);
492  Array<int> i_dofs, b_dofs, e2f;
493 
494  {
495  int info, vdim = fes->GetVDim();
496  Array<int> lvdofs;
497  Mesh *mesh = fes->GetMesh();
498  mesh->GetBdrElementAdjacentElement(bdr_el, el, info);
499  e2f.SetSize(hat_offsets[el+1]-hat_offsets[el], -1);
500  lvdofs.Reserve(A.Height());
502  mesh->Dimension()-1, info, lvdofs);
503  // Convert local element dofs to local element vdofs.
504  Ordering::DofsToVDofs<Ordering::byNODES>(e2f.Size()/vdim, vdim, lvdofs);
505  MFEM_ASSERT(lvdofs.Size() == A.Height(), "internal error");
506 #ifdef MFEM_DEBUG
507  fes->GetElementVDofs(el, vdofs);
508  for (int i = 0; i < lvdofs.Size(); i++)
509  {
510  int bd = lvdofs[i];
511  bd = (bd >= 0) ? vdofs[bd] : -1-vdofs[-1-bd];
512  MFEM_ASSERT(bvdofs[i] == bd, "internal error");
513  }
514 #endif
515  B.AdjustDofDirection(lvdofs);
517  // Create a map from local element vdofs to local boundary (face) vdofs.
518  for (int i = 0; i < lvdofs.Size(); i++)
519  {
520  e2f[lvdofs[i]] = i;
521  }
522  }
523 
524  GetIBDofs(el, i_dofs, b_dofs);
525 
526  DenseMatrix A_ii(Af_data + Af_offsets[el], i_dofs.Size(), i_dofs.Size());
527  DenseMatrix A_ib(A_ii.Data() + i_dofs.Size()*i_dofs.Size(),
528  i_dofs.Size(), b_dofs.Size());
529  DenseMatrix A_bi(A_ib.Data() + i_dofs.Size()*b_dofs.Size(),
530  b_dofs.Size(), i_dofs.Size());
531  DenseMatrix A_bb(A_bi.Data() + b_dofs.Size()*i_dofs.Size(),
532  b_dofs.Size(), b_dofs.Size());
533 
534  for (int j = 0; j < i_dofs.Size(); j++)
535  {
536  int j_f = e2f[i_dofs[j]];
537  if (j_f == -1) { continue; }
538  for (int i = 0; i < i_dofs.Size(); i++)
539  {
540  int i_f = e2f[i_dofs[i]];
541  if (i_f == -1) { continue; }
542  A_ii(i,j) += B(i_f,j_f);
543  }
544  for (int i = 0; i < b_dofs.Size(); i++)
545  {
546  int i_f = e2f[b_dofs[i]];
547  if (i_f == -1) { continue; }
548  A_bi(i,j) += B(i_f,j_f);
549  }
550  }
551  for (int j = 0; j < b_dofs.Size(); j++)
552  {
553  int j_f = e2f[b_dofs[j]];
554  if (j_f == -1) { continue; }
555  for (int i = 0; i < i_dofs.Size(); i++)
556  {
557  int i_f = e2f[i_dofs[i]];
558  if (i_f == -1) { continue; }
559  A_ib(i,j) += B(i_f,j_f);
560  }
561  for (int i = 0; i < b_dofs.Size(); i++)
562  {
563  int i_f = e2f[b_dofs[i]];
564  if (i_f == -1) { continue; }
565  A_bb(i,j) += B(i_f,j_f);
566  }
567  }
568 }
569 
571 {
572  const int skip_zeros = 1;
573  Array<int> c_dof_marker(Ct->Width());
574  Array<int> b_dofs, c_dofs;
575  const int NE = fes->GetNE();
576  DenseMatrix Cb_t, Sb_inv_Cb_t, Hb;
577 #ifndef MFEM_USE_MPI
578  H = new SparseMatrix(Ct->Width());
579 #else
580  H = pC ? NULL : new SparseMatrix(Ct->Width());
581  // V = Sb^{-1} Cb^T, for parallel non-conforming meshes
582  SparseMatrix *V = pC ? new SparseMatrix(Ct->Height(), Ct->Width()) : NULL;
583 #endif
584 
585  c_dof_marker = -1;
586  int c_mark_start = 0;
587  for (int el = 0; el < NE; el++)
588  {
589  int i_dofs_size;
590  GetBDofs(el, i_dofs_size, b_dofs);
591 
592  LUFactors LU_ii(Af_data + Af_offsets[el], Af_ipiv + Af_f_offsets[el]);
593  double *A_ib_data = LU_ii.data + i_dofs_size*i_dofs_size;
594  double *A_bi_data = A_ib_data + i_dofs_size*b_dofs.Size();
595  LUFactors LU_bb(A_bi_data + i_dofs_size*b_dofs.Size(),
596  LU_ii.ipiv + i_dofs_size);
597 
598  LU_ii.Factor(i_dofs_size);
599  LU_ii.BlockFactor(i_dofs_size, b_dofs.Size(),
600  A_ib_data, A_bi_data, LU_bb.data);
601  LU_bb.Factor(b_dofs.Size());
602 
603  // Extract Cb_t from Ct, define c_dofs
604  c_dofs.SetSize(0);
605  for (int i = 0; i < b_dofs.Size(); i++)
606  {
607  const int row = b_dofs[i];
608  const int ncols = Ct->RowSize(row);
609  const int *cols = Ct->GetRowColumns(row);
610  for (int j = 0; j < ncols; j++)
611  {
612  const int c_dof = cols[j];
613  if (c_dof_marker[c_dof] < c_mark_start)
614  {
615  c_dof_marker[c_dof] = c_mark_start + c_dofs.Size();
616  c_dofs.Append(c_dof);
617  }
618  }
619  }
620  Cb_t.SetSize(b_dofs.Size(), c_dofs.Size());
621  Cb_t = 0.0;
622  for (int i = 0; i < b_dofs.Size(); i++)
623  {
624  const int row = b_dofs[i];
625  const int ncols = Ct->RowSize(row);
626  const int *cols = Ct->GetRowColumns(row);
627  const double *vals = Ct->GetRowEntries(row);
628  for (int j = 0; j < ncols; j++)
629  {
630  const int loc_j = c_dof_marker[cols[j]] - c_mark_start;
631  Cb_t(i,loc_j) = vals[j];
632  }
633  }
634 
635  // Compute Hb = Cb Sb^{-1} Cb^t
636  Sb_inv_Cb_t = Cb_t;
637  LU_bb.Solve(Cb_t.Height(), Cb_t.Width(), Sb_inv_Cb_t.Data());
638 #ifdef MFEM_USE_MPI
639  if (!pC)
640 #endif
641  {
642  Hb.SetSize(Cb_t.Width());
643  MultAtB(Cb_t, Sb_inv_Cb_t, Hb);
644 
645  // Assemble Hb into H
646  H->AddSubMatrix(c_dofs, c_dofs, Hb, skip_zeros);
647  }
648 #ifdef MFEM_USE_MPI
649  else
650  {
651  V->AddSubMatrix(b_dofs, c_dofs, Sb_inv_Cb_t, skip_zeros);
652  }
653 #endif
654 
655  c_mark_start += c_dofs.Size();
656  MFEM_VERIFY(c_mark_start >= 0, "overflow"); // check for overflow
657  }
658  const bool fix_empty_rows = true;
659 #ifndef MFEM_USE_MPI
660  H->Finalize(skip_zeros, fix_empty_rows);
661 #else
662  ParFiniteElementSpace *c_pfes = dynamic_cast<ParFiniteElementSpace*>(c_fes);
663  if (!pC)
664  {
665  H->Finalize(skip_zeros, fix_empty_rows);
666  if (!c_pfes) { return; }
667 
668  OperatorHandle pP(pH.Type()), dH(pH.Type());
669  // TODO - construct P_pc / Dof_TrueDof_Matrix directly in the pH format
670  pP.ConvertFrom(P_pc ? P_pc : c_pfes->Dof_TrueDof_Matrix());
671  dH.MakeSquareBlockDiag(c_pfes->GetComm(),c_pfes->GlobalVSize(),
672  c_pfes->GetDofOffsets(), H);
673  pH.MakePtAP(dH, pP);
674  delete H;
675  H = NULL;
676  }
677  else
678  {
679  // TODO: add ones on the diagonal of zero rows
680  V->Finalize();
681  Array<HYPRE_BigInt> V_J(V->NumNonZeroElems());
682  MFEM_ASSERT(c_pfes, "");
683  const int c_vsize = c_fes->GetVSize();
684  HYPRE_BigInt c_ldof_offset = c_pfes->GetMyDofOffset();
685  const HYPRE_BigInt *c_face_nbr_glob_ldof = c_pfes->GetFaceNbrGlobalDofMap();
686  int *J = V->GetJ();
687  for (int i = 0; i < V_J.Size(); i++)
688  {
689  V_J[i] = J[i] < c_vsize ?
690  J[i] + c_ldof_offset :
691  c_face_nbr_glob_ldof[J[i] - c_vsize];
692  }
693  // TODO - lpH directly in the pH format
694  HypreParMatrix *lpH;
695  {
696  HypreParMatrix pV(c_pfes->GetComm(), V->Height(),
698  V->GetI(), V_J.GetData(), V->GetData(),
699  pC->ColPart(), pC->RowPart());
700  // The above constructor makes copies of all input arrays, so we can
701  // safely delete V_J and V:
702  V_J.DeleteAll();
703  delete V;
704  lpH = ParMult(pC, &pV);
705  }
706  OperatorHandle pP(pH.Type()), plpH(pH.Type());
707  // TODO - construct P_pc directly in the pH format
708  pP.ConvertFrom(P_pc);
709  plpH.ConvertFrom(lpH);
710  MFEM_VERIFY(pH.Type() != Operator::PETSC_MATIS, "To be implemented");
711  pH.MakePtAP(plpH, pP);
712  delete lpH;
713  }
714 #endif
715 }
716 
718 {
719 #ifndef MFEM_USE_MPI
720  if (!H) { ComputeH(); }
721 #else
722  if (!H && !pH.Ptr()) { ComputeH(); }
723 #endif
724 }
725 
726 void Hybridization::MultAfInv(const Vector &b, const Vector &lambda, Vector &bf,
727  int mode) const
728 {
729  // b1 = Rf^t b (assuming that Ref = 0)
730  Vector b1;
731  const SparseMatrix *R = fes->GetRestrictionMatrix();
732  if (!R)
733  {
734  b1.SetDataAndSize(b.GetData(), b.Size());
735  }
736  else
737  {
738  b1.SetSize(fes->GetVSize());
739  R->MultTranspose(b, b1);
740  }
741 
742  const int NE = fes->GetMesh()->GetNE();
743  Array<int> vdofs, i_dofs, b_dofs;
744  Vector el_vals, bf_i, i_vals, b_vals;
745  bf.SetSize(hat_offsets[NE]);
746  if (mode == 1)
747  {
748 #ifdef MFEM_USE_MPI
749  ParFiniteElementSpace *c_pfes =
750  dynamic_cast<ParFiniteElementSpace*>(c_fes);
751  if (!c_pfes)
752  {
753  Ct->Mult(lambda, bf);
754  }
755  else
756  {
757  Vector L(c_pfes->GetVSize());
758  (P_pc ? P_pc : c_pfes->GetProlongationMatrix())->Mult(lambda, L);
759  pC ? pC->MultTranspose(L, bf) : Ct->Mult(L, bf);
760  }
761 #else
762  Ct->Mult(lambda, bf);
763 #endif
764  }
765  // Apply Af^{-1}
766  Array<bool> vdof_marker(b1.Size());
767  vdof_marker = false;
768  for (int i = 0; i < NE; i++)
769  {
770  fes->GetElementVDofs(i, vdofs);
771  b1.GetSubVector(vdofs, el_vals);
772  for (int j = 0; j < vdofs.Size(); j++)
773  {
774  int vdof = vdofs[j];
775  if (vdof < 0) { vdof = -1 - vdof; }
776  if (vdof_marker[vdof]) { el_vals(j) = 0.0; }
777  else { vdof_marker[vdof] = true; }
778  }
779  bf_i.MakeRef(bf, hat_offsets[i], vdofs.Size());
780  if (mode == 1)
781  {
782  el_vals -= bf_i;
783  }
784  GetIBDofs(i, i_dofs, b_dofs);
785  el_vals.GetSubVector(i_dofs, i_vals);
786  el_vals.GetSubVector(b_dofs, b_vals);
787 
788  LUFactors LU_ii(Af_data + Af_offsets[i], Af_ipiv + Af_f_offsets[i]);
789  double *U_ib = LU_ii.data + i_dofs.Size()*i_dofs.Size();
790  double *L_bi = U_ib + i_dofs.Size()*b_dofs.Size();
791  LUFactors LU_bb(L_bi + b_dofs.Size()*i_dofs.Size(),
792  LU_ii.ipiv + i_dofs.Size());
793  LU_ii.BlockForwSolve(i_dofs.Size(), b_dofs.Size(), 1, L_bi,
794  i_vals.GetData(), b_vals.GetData());
795  LU_bb.Solve(b_dofs.Size(), 1, b_vals.GetData());
796  bf_i = 0.0;
797  if (mode == 1)
798  {
799  LU_ii.BlockBackSolve(i_dofs.Size(), b_dofs.Size(), 1, U_ib,
800  b_vals.GetData(), i_vals.GetData());
801  bf_i.SetSubVector(i_dofs, i_vals);
802  }
803  bf_i.SetSubVector(b_dofs, b_vals);
804  }
805 }
806 
807 void Hybridization::ReduceRHS(const Vector &b, Vector &b_r) const
808 {
809  // bf = Af^{-1} Rf^t b
810  Vector bf;
811  MultAfInv(b, b, bf, 0);
812 
813  // b_r = Cf bf
814 #ifdef MFEM_USE_MPI
815  ParFiniteElementSpace *c_pfes = dynamic_cast<ParFiniteElementSpace*>(c_fes);
816  if (!c_pfes)
817  {
818  b_r.SetSize(Ct->Width());
819  Ct->MultTranspose(bf, b_r);
820  }
821  else
822  {
823  Vector bl(pC ? pC->Height() : Ct->Width());
824  if (pC)
825  {
826  pC->Mult(bf, bl);
827  }
828  else
829  {
830  Ct->MultTranspose(bf, bl);
831  }
832  b_r.SetSize(pH.Ptr()->Height());
833  (P_pc ? P_pc : c_pfes->GetProlongationMatrix())->MultTranspose(bl, b_r);
834  }
835 #else
836  b_r.SetSize(Ct->Width());
837  Ct->MultTranspose(bf, b_r);
838 #endif
839 }
840 
841 void Hybridization::ComputeSolution(const Vector &b, const Vector &sol_r,
842  Vector &sol) const
843 {
844  // bf = Af^{-1} ( Rf^t - Cf^t sol_r )
845  Vector bf;
846  MultAfInv(b, sol_r, bf, 1);
847 
848  // sol = Rf bf
849  GridFunction s;
850  const SparseMatrix *R = fes->GetRestrictionMatrix();
851  if (!R)
852  {
853  MFEM_ASSERT(sol.Size() == fes->GetVSize(), "");
854  s.MakeRef(fes, sol, 0);
855  }
856  else
857  {
858  s.SetSpace(fes);
859  R->MultTranspose(sol, s);
860  }
861  const int NE = fes->GetMesh()->GetNE();
862  Array<int> vdofs;
863  for (int i = 0; i < NE; i++)
864  {
865  fes->GetElementVDofs(i, vdofs);
866  for (int j = hat_offsets[i]; j < hat_offsets[i+1]; j++)
867  {
868  if (hat_dofs_marker[j] == 1) { continue; } // skip essential b.c.
869  int vdof = vdofs[j-hat_offsets[i]];
870  if (vdof >= 0) { s(vdof) = bf(j); }
871  else { s(-1-vdof) = -bf(j); }
872  }
873  }
874  if (R)
875  {
876  R->Mult(s, sol); // assuming that Ref = 0
877  }
878 }
879 
881 {
882  delete H;
883  H = NULL;
884 #ifdef MFEM_USE_MPI
885  pH.Clear();
886 #endif
887 }
888 
889 }
Abstract class for all finite elements.
Definition: fe_base.hpp:235
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
Array< int > Af_offsets
int RowSize(const int i) const
Returns the number of elements in row i.
Definition: sparsemat.cpp:344
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
void ConvertFrom(OperatorHandle &A)
Convert the given OperatorHandle A to the currently set type id.
Definition: handle.cpp:200
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
const SparseMatrix * GetConformingProlongation() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition: fespace.cpp:1232
FiniteElementSpace * c_fes
Array< int > Af_f_offsets
ParMesh * GetParMesh() const
Definition: pfespace.hpp:277
Array< Slave > slaves
Definition: ncmesh.hpp:231
virtual void Finalize(int skip_zeros=1)
Finalize the matrix initialization, switching the storage format from LIL to CSR. ...
Definition: sparsemat.hpp:545
virtual const Operator * GetProlongationMatrix() const
The returned Operator is owned by the FiniteElementSpace.
Definition: pfespace.cpp:1153
void GetBDofs(int el, int &num_idofs, Array< int > &b_dofs) const
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:513
HypreParMatrix * P_pc
int TrueVSize() const
Obsolete, kept for backward compatibility.
Definition: pfespace.hpp:434
int * GetRowColumns(const int row)
Return a pointer to the column indices in a row.
Definition: sparsemat.cpp:391
void BooleanMult(const Array< int > &x, Array< int > &y) const
y = A * x, treating all entries as booleans (zero=false, nonzero=true).
Definition: sparsemat.cpp:1022
Lists all edges/faces in the nonconforming mesh.
Definition: ncmesh.hpp:227
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
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
Operator::Type Type() const
Get the currently set operator type id.
Definition: handle.hpp:99
void BlockFactor(int m, int n, double *A12, double *A21, double *A22) const
Definition: densemat.cpp:3530
void BlockBackSolve(int m, int n, int r, const double *U12, const double *X2, double *Y1) const
Definition: densemat.cpp:3565
A specialized ElementTransformation class representing a face and its two neighboring elements...
Definition: eltrans.hpp:480
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition: vector.cpp:547
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
int * GetI()
Return the array I.
Definition: sparsemat.hpp:222
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:923
Abstract parallel finite element space.
Definition: pfespace.hpp:28
const FiniteElement * GetFaceElement(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th face in the ...
Definition: fespace.cpp:3133
HypreParMatrix * GetPartialConformingInterpolation()
For a non-conforming mesh, construct and return the interpolation matrix from the partially conformin...
Definition: pfespace.cpp:964
FaceElementTransformations * GetInteriorFaceTransformations(int FaceNo)
Definition: mesh.hpp:1268
double * GetRowEntries(const int row)
Return a pointer to the entries in a row.
Definition: sparsemat.cpp:405
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:209
HypreParMatrix * ParMult(const HypreParMatrix *A, const HypreParMatrix *B, bool own_matrix)
Definition: hypre.cpp:2849
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:1069
int GetConformingVSize() const
Definition: fespace.hpp:596
void Reset()
Destroy the current hybridization matrix while preserving the computed constraint matrix and the set ...
void DeleteAll()
Delete the whole array.
Definition: array.hpp:846
const NCList & GetSharedList(int entity)
Helper to get shared vertices/edges/faces (&#39;entity&#39; == 0/1/2 resp.).
Definition: pncmesh.hpp:128
SparseMatrix * Ct
DofTransformation * GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i&#39;th boundary element.
Definition: fespace.cpp:297
ParNCMesh * pncmesh
Definition: pmesh.hpp:388
double * GetData()
Return the element data, i.e. the array A.
Definition: sparsemat.hpp:232
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
Definition: mesh.cpp:5376
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:614
void Finalize()
Finalize the construction of the hybridized matrix.
bool Nonconforming() const
Definition: mesh.hpp:1651
HYPRE_BigInt GlobalVSize() const
Definition: pfespace.hpp:283
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 AssembleBdrMatrix(int bdr_el, const DenseMatrix &A)
Assemble the boundary element matrix A into the hybridized system matrix.
MPI_Comm GetComm() const
Definition: pfespace.hpp:273
Data type sparse matrix.
Definition: sparsemat.hpp:50
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:67
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
int GetSharedFace(int sface) const
Return the local face index for the given shared face.
Definition: pmesh.cpp:3146
double b
Definition: lissajous.cpp:42
FaceElementTransformations * GetFaceElementTransformations(int FaceNo, int mask=31) override
Definition: pmesh.cpp:2995
const HYPRE_BigInt * GetFaceNbrGlobalDofMap()
Definition: pfespace.hpp:399
void ComputeSolution(const Vector &b, const Vector &sol_r, Vector &sol) const
void BlockForwSolve(int m, int n, int r, const double *L21, double *B1, double *B2) const
Definition: densemat.cpp:3556
Array< int > hat_offsets
FaceElementTransformations * GetSharedFaceTransformations(int sf, bool fill2=true)
Definition: pmesh.cpp:3011
HypreParMatrix * Dof_TrueDof_Matrix() const
The true dof-to-dof interpolation matrix.
Definition: pfespace.hpp:321
void Reserve(int capacity)
Ensures that the allocated size is at least the given size.
Definition: array.hpp:156
void SubDofOrder(Geometry::Type Geom, int SDim, int Info, Array< int > &dofs) const
Get the local dofs for a given sub-manifold.
Definition: fe_coll.cpp:397
void GetBdrElementAdjacentElement(int bdr_el, int &el, int &info) const
For the given boundary element, bdr_el, return its adjacent element and its info, i...
Definition: mesh.cpp:6206
HypreParMatrix * Transpose() const
Returns the transpose of *this.
Definition: hypre.cpp:1607
virtual void MakeRef(FiniteElementSpace *f, double *v)
Make the GridFunction reference external data on a new FiniteElementSpace.
Definition: gridfunc.cpp:204
HYPRE_BigInt GetMyDofOffset() const
Definition: pfespace.cpp:1143
double * data
Definition: densemat.hpp:599
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:581
virtual void PrintMatlab(std::ostream &out=mfem::out) const
Prints matrix in matlab format.
Definition: sparsemat.cpp:3288
HYPRE_BigInt * GetDofOffsets() const
Definition: pfespace.hpp:281
int Dimension() const
Definition: mesh.hpp:1006
void MultAfInv(const Vector &b, const Vector &lambda, Vector &bf, int mode) const
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
HYPRE_BigInt GetGlobalNumCols() const
Return the global number of columns.
Definition: hypre.hpp:639
void Threshold(double eps)
Replace small entries, abs(a_ij) &lt;= eps, with zero.
Definition: densemat.cpp:2169
BilinearFormIntegrator * c_bfi
void BooleanMult(int alpha, const int *x, int beta, int *y)
The &quot;Boolean&quot; analog of y = alpha * A * x + beta * y, where elements in the sparsity pattern of the m...
Definition: hypre.hpp:714
double MaxMaxNorm() const
Compute the norm ||A|| = max_{ij} |A_{ij}|.
Definition: densemat.cpp:808
double * Data() const
Returns the matrix data array.
Definition: densemat.hpp:111
SparseMatrix * H
ComplexDenseMatrix * MultAtB(const ComplexDenseMatrix &A, const ComplexDenseMatrix &B)
Multiply the complex conjugate transpose of a matrix A with a matrix B. A^H*B.
void AssembleMatrix(int el, const DenseMatrix &A)
Assemble the element matrix A into the hybridized system matrix.
MPI_Comm GetComm() const
Definition: pmesh.hpp:351
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Definition: sparsemat.cpp:729
void GetIBDofs(int el, Array< int > &i_dofs, Array< int > &b_dofs) const
MFEM_HOST_DEVICE void MultTranspose(const int height, const int width, const TA *data, const TX *x, TY *y)
Matrix transpose vector multiplication: y = At x, where the matrix A is of size height x width with g...
Definition: kernels.hpp:196
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:679
HYPRE_Int HYPRE_BigInt
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
Definition: table.cpp:413
OperatorHandle pH
void Clear()
Clear the OperatorHandle, deleting the held Operator (if owned), while leaving the type id unchanged...
Definition: handle.hpp:124
~Hybridization()
Destructor.
void SetDataAndSize(double *d, int s)
Set the Vector data and size.
Definition: vector.hpp:157
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 AssembleFaceMatrix(const FiniteElement &el1, const FiniteElement &el2, FaceElementTransformations &Trans, DenseMatrix &elmat)
Definition: bilininteg.cpp:139
HYPRE_BigInt GetGlobalNumRows() const
Return the global number of rows.
Definition: hypre.hpp:635
HYPRE_BigInt * ColPart()
Returns the column partitioning.
Definition: hypre.hpp:573
FiniteElementSpace * fes
HYPRE_BigInt * RowPart()
Returns the row partitioning.
Definition: hypre.hpp:569
HypreParMatrix * pC
int dim
Definition: ex24.cpp:53
virtual bool Factor(int m, double TOL=0.0)
Compute the LU factorization of the current matrix.
Definition: densemat.cpp:3230
void MultTranspose(const Vector &x, Vector &y) const
Multiply a vector with the transposed matrix. y = At * x.
Definition: sparsemat.cpp:904
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th element in t...
Definition: fespace.cpp:2781
const FiniteElement * GetFaceNbrFaceFE(int i) const
Definition: pfespace.cpp:1531
virtual void SetSpace(FiniteElementSpace *f)
Associate a new FiniteElementSpace with the GridFunction.
Definition: gridfunc.cpp:196
int GetFaceNbrVSize() const
Definition: pfespace.hpp:394
Array< int > hat_dofs_marker
void GetFaceVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i&#39;th face element (2D and 3D).
Definition: fespace.cpp:312
T & Last()
Return the last element in the array.
Definition: array.hpp:784
void GetFaceNbrFaceVDofs(int i, Array< int > &vdofs) const
Definition: pfespace.cpp:1489
void MakeRef(T *, int)
Make this Array a reference to a pointer.
Definition: array.hpp:864
Vector data type.
Definition: vector.hpp:60
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:601
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
Definition: vector.hpp:577
ID for class HypreParMatrix.
Definition: operator.hpp:260
HYPRE_Int Mult(HypreParVector &x, HypreParVector &y, double alpha=1.0, double beta=0.0) const
Computes y = alpha * A * x + beta * y.
Definition: hypre.cpp:1732
void ReduceRHS(const Vector &b, Vector &b_r) const
RefCoord s[3]
void Init(const Array< int > &ess_tdof_list)
Prepare the Hybridization object for assembly.
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:66
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
virtual const SparseMatrix * GetRestrictionMatrix() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition: fespace.hpp:495
Hybridization(FiniteElementSpace *fespace, FiniteElementSpace *c_fespace)
Constructor.
Class for parallel meshes.
Definition: pmesh.hpp:32
void AdjustDofDirection(Array< int > &dofs)
Definition: densemat.cpp:2096
ID for class PetscParMatrix, MATIS format.
Definition: operator.hpp:262
void GenerateOffsets(int N, HYPRE_BigInt loc_sizes[], Array< HYPRE_BigInt > *offsets[]) const
Definition: pmesh.cpp:1933
void MakePtAP(OperatorHandle &A, OperatorHandle &P)
Reset the OperatorHandle to hold the product P^t A P.
Definition: handle.cpp:123
static void AdjustVDofs(Array< int > &vdofs)
Definition: fespace.cpp:267
void SetType(Operator::Type tid)
Invoke Clear() and set a new type id.
Definition: handle.hpp:132