MFEM  v3.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, 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 "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),
33  Af_data(NULL), Af_ipiv(NULL)
34 {
35 #ifdef MFEM_USE_MPI
36  pH = NULL;
37 #endif
38 }
39 
41 {
42 #ifdef MFEM_USE_MPI
43  delete pH;
44 #endif
45  delete [] Af_ipiv;
46  delete [] Af_data;
47  delete H;
48  delete Ct;
49  delete c_bfi;
50 }
51 
52 void Hybridization::Init(const Array<int> &ess_tdof_list)
53 {
54  if (Ct) { return; }
55 
56  // Assemble the constraint matrix C
57  // TODO: ParNCMesh
58 
59  // count the number of dofs in the discontinuous version of fes:
60  const int NE = fes->GetNE();
61  Array<int> vdofs, c_vdofs;
62  int num_hat_dofs = 0;
63  hat_offsets.SetSize(NE+1);
64  hat_offsets[0] = 0;
65  for (int i = 0; i < NE; i++)
66  {
67  fes->GetElementVDofs(i, vdofs);
68  num_hat_dofs += vdofs.Size();
69  hat_offsets[i+1] = num_hat_dofs;
70  }
71 
72 #ifdef MFEM_USE_MPI
73  ParFiniteElementSpace *c_pfes =
74  dynamic_cast<ParFiniteElementSpace*>(c_fes);
75  ParMesh *pmesh = c_pfes ? c_pfes->GetParMesh() : NULL;
76  MFEM_VERIFY(!c_pfes || c_pfes->Conforming(),
77  "parallel non-conforming AMR meshes are not supported yet");
78 #endif
79 
80  Ct = new SparseMatrix(num_hat_dofs, c_fes->GetVSize());
81 
82  if (c_bfi)
83  {
84  const int skip_zeros = 1;
85  DenseMatrix elmat;
87  Mesh *mesh = fes->GetMesh();
88  int num_faces = mesh->GetNumFaces();
89  for (int i = 0; i < num_faces; i++)
90  {
91  FTr = mesh->GetInteriorFaceTransformations(i);
92  if (!FTr) { continue; }
93 
94  int o1 = hat_offsets[FTr->Elem1No];
95  int s1 = hat_offsets[FTr->Elem1No+1] - o1;
96  int o2 = hat_offsets[FTr->Elem2No];
97  int s2 = hat_offsets[FTr->Elem2No+1] - o2;
98  vdofs.SetSize(s1 + s2);
99  for (int j = 0; j < s1; j++)
100  {
101  vdofs[j] = o1 + j;
102  }
103  for (int j = 0; j < s2; j++)
104  {
105  vdofs[s1+j] = o2 + j;
106  }
107  c_fes->GetFaceVDofs(i, c_vdofs);
109  *fes->GetFE(FTr->Elem1No),
110  *fes->GetFE(FTr->Elem2No),
111  *FTr, elmat);
112  // zero-out small elements in elmat
113  elmat.Threshold(1e-12 * elmat.MaxMaxNorm());
114  Ct->AddSubMatrix(vdofs, c_vdofs, elmat, skip_zeros);
115  }
116 #ifdef MFEM_USE_MPI
117  if (pmesh)
118  {
119  // Assemble local contribution to Ct from shared faces
120  const int num_shared_faces = pmesh->GetNSharedFaces();
121  for (int i = 0; i < num_shared_faces; i++)
122  {
123  int face_no = pmesh->GetSharedFace(i);
124  FTr = pmesh->GetFaceElementTransformations(face_no);
125  MFEM_ASSERT(FTr->Elem2No < 0, "");
126 
127  int o1 = hat_offsets[FTr->Elem1No];
128  int s1 = hat_offsets[FTr->Elem1No+1] - o1;
129  vdofs.SetSize(s1);
130  for (int j = 0; j < s1; j++)
131  {
132  vdofs[j] = o1 + j;
133  }
134  c_fes->GetFaceVDofs(face_no, c_vdofs);
135  const FiniteElement *fe = fes->GetFE(FTr->Elem1No);
137  *fe, *fe, *FTr, elmat);
138  // zero-out small elements in elmat
139  elmat.Threshold(1e-12 * elmat.MaxMaxNorm());
140  Ct->AddSubMatrix(vdofs, c_vdofs, elmat, skip_zeros);
141  }
142  }
143 #endif
144  Ct->Finalize(skip_zeros);
145  }
146  else
147  {
148  // Check if c_fes is really needed here.
149  MFEM_ABORT("TODO: algebraic definition of C");
150  }
151 
152 #ifdef MFEM_DEBUG_HYBRIDIZATION_CP
153  // Debug: write C and P to file
154  {
155  std::ofstream C_file("C_matrix.txt");
156  SparseMatrix *C = Transpose(*Ct);
157  C->PrintMatlab(C_file);
158  delete C;
159 
161  if (P)
162  {
163  std::ofstream P_file("P_matrix.txt");
164  P->PrintMatlab(P_file);
165  }
166  }
167 #endif
168 
169  // Define the "free" (0) and "essential" (1) hat_dofs.
170  // The "essential" hat_dofs are those that depend only on essential cdofs;
171  // all other hat_dofs are "free".
172  hat_dofs_marker.SetSize(num_hat_dofs);
173  Array<int> free_tdof_marker;
174 #ifdef MFEM_USE_MPI
175  ParFiniteElementSpace *pfes = dynamic_cast<ParFiniteElementSpace*>(fes);
176  free_tdof_marker.SetSize(pfes ? pfes->TrueVSize() :
178 #else
179  free_tdof_marker.SetSize(fes->GetConformingVSize());
180 #endif
181  free_tdof_marker = 1;
182  for (int i = 0; i < ess_tdof_list.Size(); i++)
183  {
184  free_tdof_marker[ess_tdof_list[i]] = 0;
185  }
186  Array<int> free_vdofs_marker;
187 #ifdef MFEM_USE_MPI
188  if (!pfes)
189  {
191  if (!cP)
192  {
193  free_vdofs_marker.MakeRef(free_tdof_marker);
194  }
195  else
196  {
197  free_vdofs_marker.SetSize(fes->GetVSize());
198  cP->BooleanMult(free_tdof_marker, free_vdofs_marker);
199  }
200  }
201  else
202  {
203  HypreParMatrix *P = pfes->Dof_TrueDof_Matrix();
204  free_vdofs_marker.SetSize(fes->GetVSize());
205  P->BooleanMult(1, free_tdof_marker, 0, free_vdofs_marker);
206  }
207 #else
209  if (!cP)
210  {
211  free_vdofs_marker.MakeRef(free_tdof_marker);
212  }
213  else
214  {
215  free_vdofs_marker.SetSize(fes->GetVSize());
216  cP->BooleanMult(free_tdof_marker, free_vdofs_marker);
217  }
218 #endif
219  for (int i = 0; i < NE; i++)
220  {
221  fes->GetElementVDofs(i, vdofs);
223  for (int j = 0; j < vdofs.Size(); j++)
224  {
225  hat_dofs_marker[hat_offsets[i]+j] = ! free_vdofs_marker[vdofs[j]];
226  }
227  }
228 #ifndef MFEM_DEBUG
229  // In DEBUG mode this array is used below.
230  free_tdof_marker.DeleteAll();
231 #endif
232  free_vdofs_marker.DeleteAll();
233  // Split the "free" (0) hat_dofs into "internal" (0) or "boundary" (-1).
234  // The "internal" hat_dofs are those "free" hat_dofs for which the
235  // corresponding column in C is zero; otherwise the free hat_dof is
236  // "boundary".
237  for (int i = 0; i < num_hat_dofs; i++)
238  {
239  // skip "essential" hat_dofs and empty rows in Ct
240  if (hat_dofs_marker[i] != 1 && Ct->RowSize(i) > 0)
241  {
242  hat_dofs_marker[i] = -1; // mark this hat_dof as "boundary"
243  }
244  }
245 
246  H = new SparseMatrix(Ct->Width());
247 
248  // Define Af_offsets and Af_f_offsets
249  Af_offsets.SetSize(NE+1);
250  Af_offsets[0] = 0;
251  Af_f_offsets.SetSize(NE+1);
252  Af_f_offsets[0] = 0;
253  // #define MFEM_DEBUG_HERE // uncomment to enable printing of hat dofs stats
254 #ifdef MFEM_DEBUG_HERE
255  int b_size = 0;
256 #endif
257  for (int i = 0; i < NE; i++)
258  {
259  int f_size = 0; // count the "free" hat_dofs in element i
260  for (int j = hat_offsets[i]; j < hat_offsets[i+1]; j++)
261  {
262  if (hat_dofs_marker[j] != 1) { f_size++; }
263 #ifdef MFEM_DEBUG_HERE
264  if (hat_dofs_marker[j] == -1) { b_size++; }
265 #endif
266  }
267  Af_offsets[i+1] = Af_offsets[i] + f_size*f_size;
268  Af_f_offsets[i+1] = Af_f_offsets[i] + f_size;
269  }
270 
271 #ifdef MFEM_DEBUG_HERE
272 #ifndef MFEM_USE_MPI
273  int myid = 0;
274 #else
275  int myid = pmesh ? pmesh->GetMyRank() : 0;
276 #endif
277  int i_size = Af_f_offsets[NE] - b_size;
278  int e_size = num_hat_dofs - (i_size + b_size);
279  std::cout << "\nHybridization::Init:"
280  << " [" << myid << "] hat dofs - \"internal\": " << i_size
281  << ", \"boundary\": " << b_size
282  << ", \"essential\": " << e_size << '\n' << std::endl;
283 #undef MFEM_DEBUG_HERE
284 #endif
285 
286  Af_data = new double[Af_offsets[NE]];
287  Af_ipiv = new int[Af_f_offsets[NE]];
288 
289 #ifdef MFEM_DEBUG
290  // check that Ref = 0
291  const SparseMatrix *R = fes->GetRestrictionMatrix();
292  if (!R) { return; }
293  Array<int> vdof_marker(fes->GetVSize()); // 0 - f, 1 - e
294  vdof_marker = 0;
295  for (int i = 0; i < NE; i++)
296  {
297  fes->GetElementVDofs(i, vdofs);
299  for (int j = 0; j < vdofs.Size(); j++)
300  {
301  if (hat_dofs_marker[hat_offsets[i]+j] == 1) // "essential" hat dof
302  {
303  vdof_marker[vdofs[j]] = 1;
304  }
305  }
306  }
307  for (int tdof = 0; tdof < R->Height(); tdof++)
308  {
309  if (free_tdof_marker[tdof]) { continue; }
310 
311  const int ncols = R->RowSize(tdof);
312  const int *cols = R->GetRowColumns(tdof);
313  const double *vals = R->GetRowEntries(tdof);
314  for (int j = 0; j < ncols; j++)
315  {
316  if (std::abs(vals[j]) != 0.0 && vdof_marker[cols[j]] == 0)
317  {
318  MFEM_ABORT("Ref != 0");
319  }
320  }
321  }
322 #endif
323 }
324 
326  int el, Array<int> &i_dofs, Array<int> &b_dofs) const
327 {
328  int h_start, h_end;
329 
330  h_start = hat_offsets[el];
331  h_end = hat_offsets[el+1];
332  i_dofs.Reserve(h_end-h_start);
333  i_dofs.SetSize(0);
334  b_dofs.Reserve(h_end-h_start);
335  b_dofs.SetSize(0);
336  for (int i = h_start; i < h_end; i++)
337  {
338  int mark = hat_dofs_marker[i];
339  if (mark == 0) { i_dofs.Append(i-h_start); }
340  else if (mark == -1) { b_dofs.Append(i-h_start); }
341  }
342 }
343 
345 {
346  Array<int> i_dofs, b_dofs;
347 
348  GetIBDofs(el, i_dofs, b_dofs);
349 
350  DenseMatrix A_ii(Af_data + Af_offsets[el], i_dofs.Size(), i_dofs.Size());
351  DenseMatrix A_ib(A_ii.Data() + i_dofs.Size()*i_dofs.Size(),
352  i_dofs.Size(), b_dofs.Size());
353  DenseMatrix A_bi(A_ib.Data() + i_dofs.Size()*b_dofs.Size(),
354  b_dofs.Size(), i_dofs.Size());
355  DenseMatrix A_bb(A_bi.Data() + b_dofs.Size()*i_dofs.Size(),
356  b_dofs.Size(), b_dofs.Size());
357 
358  for (int j = 0; j < i_dofs.Size(); j++)
359  {
360  int j_dof = i_dofs[j];
361  for (int i = 0; i < i_dofs.Size(); i++)
362  {
363  A_ii(i,j) = A(i_dofs[i],j_dof);
364  }
365  for (int i = 0; i < b_dofs.Size(); i++)
366  {
367  A_bi(i,j) = A(b_dofs[i],j_dof);
368  }
369  }
370  for (int j = 0; j < b_dofs.Size(); j++)
371  {
372  int j_dof = b_dofs[j];
373  for (int i = 0; i < i_dofs.Size(); i++)
374  {
375  A_ib(i,j) = A(i_dofs[i],j_dof);
376  }
377  for (int i = 0; i < b_dofs.Size(); i++)
378  {
379  A_bb(i,j) = A(b_dofs[i],j_dof);
380  }
381  }
382 
383  LUFactors LU_ii(A_ii.Data(), Af_ipiv + Af_f_offsets[el]);
384  LUFactors LU_bb(A_bb.Data(), LU_ii.ipiv + i_dofs.Size());
385 
386  LU_ii.Factor(i_dofs.Size());
387  LU_ii.BlockFactor(i_dofs.Size(), b_dofs.Size(),
388  A_ib.Data(), A_bi.Data(), A_bb.Data());
389  LU_bb.Factor(b_dofs.Size());
390 
391  // Extract Cb_t from Ct
392  std::map<int,int> Cb_g2l;
393  for (int i = 0; i < b_dofs.Size(); i++)
394  {
395  const int row = hat_offsets[el] + b_dofs[i];
396  const int ncols = Ct->RowSize(row);
397  const int *cols = Ct->GetRowColumns(row);
398  for (int j = 0; j < ncols; j++)
399  {
400  const std::pair<const int,int> p(cols[j], (int)Cb_g2l.size());
401  Cb_g2l.insert(p);
402  }
403  }
404  Array<int> c_dofs((int)Cb_g2l.size());
405  for (std::map<int, int>::iterator it = Cb_g2l.begin();
406  it != Cb_g2l.end(); ++it)
407  {
408  c_dofs[it->second] = it->first;
409  }
410  DenseMatrix Cb_t(b_dofs.Size(), c_dofs.Size()); // Cb_t is init with 0
411  for (int i = 0; i < b_dofs.Size(); i++)
412  {
413  const int row = hat_offsets[el] + b_dofs[i];
414  const int ncols = Ct->RowSize(row);
415  const int *cols = Ct->GetRowColumns(row);
416  const double *vals = Ct->GetRowEntries(row);
417  for (int j = 0; j < ncols; j++)
418  {
419  Cb_t(i,Cb_g2l[cols[j]]) = vals[j];
420  }
421  }
422 
423  // Compute Hb = Cb Sb^{-1} Cb^t
424  DenseMatrix Sb_inv_Cb_t(Cb_t);
425  LU_bb.Solve(Cb_t.Height(), Cb_t.Width(), Sb_inv_Cb_t.Data());
426  DenseMatrix Hb(Cb_t.Width());
427  MultAtB(Cb_t, Sb_inv_Cb_t, Hb);
428 
429  // Assemble Hb into H
430  const int skip_zeros = 1;
431  H->AddSubMatrix(c_dofs, c_dofs, Hb, skip_zeros);
432 }
433 
435 {
436 #ifndef MFEM_USE_MPI
437  H->Finalize();
438 #else
439  if (!H) { return; }
440  H->Finalize();
441  ParFiniteElementSpace *c_pfes =
442  dynamic_cast<ParFiniteElementSpace*>(c_fes);
443  if (!c_pfes) { return; }
444  HypreParMatrix *dH =
445  new HypreParMatrix(c_pfes->GetComm(), c_pfes->GlobalVSize(),
446  c_pfes->GetDofOffsets(), H);
447  // TODO: ParNCMesh
448  pH = RAP(dH, c_pfes->Dof_TrueDof_Matrix());
449  delete dH;
450  delete H;
451  H = NULL;
452 #endif
453  // TODO: add ones on the diagonal of zero rows.
454 }
455 
456 void Hybridization::MultAfInv(const Vector &b, const Vector &lambda, Vector &bf,
457  int mode) const
458 {
459  // b1 = Rf^t b (assuming that Ref = 0)
460  Vector b1;
461  const SparseMatrix *R = fes->GetRestrictionMatrix();
462  if (!R)
463  {
464  b1.SetDataAndSize(b.GetData(), b.Size());
465  }
466  else
467  {
468  b1.SetSize(fes->GetVSize());
469  R->MultTranspose(b, b1);
470  }
471 
472  const int NE = fes->GetMesh()->GetNE();
473  Array<int> vdofs, i_dofs, b_dofs;
474  Vector el_vals, bf_i, i_vals, b_vals;
475  bf.SetSize(hat_offsets[NE]);
476  if (mode == 1)
477  {
478 #ifdef MFEM_USE_MPI
479  ParFiniteElementSpace *c_pfes =
480  dynamic_cast<ParFiniteElementSpace*>(c_fes);
481  if (!c_pfes)
482  {
483  Ct->Mult(lambda, bf);
484  }
485  else
486  {
487  // TODO: ParNCMesh
488  Vector L(c_pfes->GetVSize());
489  c_pfes->Dof_TrueDof_Matrix()->Mult(lambda, L);
490  Ct->Mult(L, bf);
491  }
492 #else
493  Ct->Mult(lambda, bf);
494 #endif
495  }
496  // Apply Af^{-1}
497  Array<bool> vdof_marker(b1.Size());
498  vdof_marker = false;
499  for (int i = 0; i < NE; i++)
500  {
501  fes->GetElementVDofs(i, vdofs);
502  b1.GetSubVector(vdofs, el_vals);
503  for (int j = 0; j < vdofs.Size(); j++)
504  {
505  int vdof = vdofs[j];
506  if (vdof < 0) { vdof = -1 - vdof; }
507  if (vdof_marker[vdof]) { el_vals(j) = 0.0; }
508  else { vdof_marker[vdof] = true; }
509  }
510  bf_i.SetDataAndSize(&bf[hat_offsets[i]], vdofs.Size());
511  if (mode == 1)
512  {
513  el_vals -= bf_i;
514  }
515  GetIBDofs(i, i_dofs, b_dofs);
516  el_vals.GetSubVector(i_dofs, i_vals);
517  el_vals.GetSubVector(b_dofs, b_vals);
518 
519  LUFactors LU_ii(Af_data + Af_offsets[i], Af_ipiv + Af_f_offsets[i]);
520  double *U_ib = LU_ii.data + i_dofs.Size()*i_dofs.Size();
521  double *L_bi = U_ib + i_dofs.Size()*b_dofs.Size();
522  LUFactors LU_bb(L_bi + b_dofs.Size()*i_dofs.Size(),
523  LU_ii.ipiv + i_dofs.Size());
524  LU_ii.BlockForwSolve(i_dofs.Size(), b_dofs.Size(), 1, L_bi,
525  i_vals.GetData(), b_vals.GetData());
526  LU_bb.Solve(b_dofs.Size(), 1, b_vals.GetData());
527  bf_i = 0.0;
528  if (mode == 1)
529  {
530  LU_ii.BlockBackSolve(i_dofs.Size(), b_dofs.Size(), 1, U_ib,
531  b_vals.GetData(), i_vals.GetData());
532  bf_i.SetSubVector(i_dofs, i_vals);
533  }
534  bf_i.SetSubVector(b_dofs, b_vals);
535  }
536 }
537 
538 void Hybridization::ReduceRHS(const Vector &b, Vector &b_r) const
539 {
540  // bf = Af^{-1} Rf^t b
541  Vector bf;
542  MultAfInv(b, b, bf, 0);
543 
544  // b_r = Cf bf
545 #ifdef MFEM_USE_MPI
546  ParFiniteElementSpace *c_pfes =
547  dynamic_cast<ParFiniteElementSpace*>(c_fes);
548  if (!c_pfes)
549  {
550  b_r.SetSize(Ct->Width());
551  Ct->MultTranspose(bf, b_r);
552  }
553  else
554  {
555  // TODO: ParNCMesh
556  Vector bl(Ct->Width());
557  Ct->MultTranspose(bf, bl);
558  b_r.SetSize(pH->Height());
559  c_pfes->Dof_TrueDof_Matrix()->MultTranspose(bl, b_r);
560  }
561 #else
562  b_r.SetSize(Ct->Width());
563  Ct->MultTranspose(bf, b_r);
564 #endif
565 }
566 
567 void Hybridization::ComputeSolution(const Vector &b, const Vector &sol_r,
568  Vector &sol) const
569 {
570  // bf = Af^{-1} ( Rf^t - Cf^t sol_r )
571  Vector bf;
572  MultAfInv(b, sol_r, bf, 1);
573 
574  // sol = Rf bf
575  GridFunction s;
576  const SparseMatrix *R = fes->GetRestrictionMatrix();
577  if (!R)
578  {
579  MFEM_ASSERT(sol.Size() == fes->GetVSize(), "");
580  s.Update(fes, sol, 0);
581  }
582  else
583  {
584  s.Update(fes);
585  R->MultTranspose(sol, s);
586  }
587  const int NE = fes->GetMesh()->GetNE();
588  Array<int> vdofs;
589  for (int i = 0; i < NE; i++)
590  {
591  fes->GetElementVDofs(i, vdofs);
592  for (int j = hat_offsets[i]; j < hat_offsets[i+1]; j++)
593  {
594  if (hat_dofs_marker[j] == 1) { continue; } // skip essential b.c.
595  int vdof = vdofs[j-hat_offsets[i]];
596  if (vdof >= 0) { s(vdof) = bf(j); }
597  else { s(-1-vdof) = -bf(j); }
598  }
599  }
600  if (R)
601  {
602  R->Mult(s, sol); // assuming that Ref = 0
603  }
604 }
605 
606 }
Abstract class for Finite Elements.
Definition: fe.hpp:44
int Size() const
Logical size of the array.
Definition: array.hpp:109
Array< int > Af_offsets
int RowSize(const int i) const
Returns the number of elements in row i.
Definition: sparsemat.cpp:190
int GetVSize() const
Definition: fespace.hpp:164
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:27
FiniteElementSpace * c_fes
Array< int > Af_f_offsets
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:963
void SetSize(int s)
Resizes the vector if the new size is different.
Definition: vector.hpp:259
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:161
int * GetRowColumns(const int row)
Return a pointer to the column indices in a row.
Definition: sparsemat.cpp:237
void BooleanMult(const Array< int > &x, Array< int > &y) const
y = A * x, but treat all elements as booleans (zero=false, nonzero=true).
Definition: sparsemat.cpp:573
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols.
Definition: operator.hpp:41
int GetNSharedFaces() const
Return the number of shared faces (3D), edges (2D), vertices (1D)
Definition: pmesh.cpp:1409
void BlockBackSolve(int m, int n, int r, const double *U12, const double *X2, double *Y1) const
Definition: densemat.cpp:3821
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Definition: vector.cpp:457
HYPRE_Int * GetDofOffsets()
Definition: pfespace.hpp:151
Data type dense matrix using column-major storage.
Definition: densemat.hpp:22
int Size() const
Returns the size of the vector.
Definition: vector.hpp:84
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:454
Abstract parallel finite element space.
Definition: pfespace.hpp:28
const FiniteElement * GetFaceElement(int i) const
Definition: fespace.cpp:1445
void Factor(int m)
Definition: densemat.cpp:3557
FaceElementTransformations * GetInteriorFaceTransformations(int FaceNo)
Definition: mesh.hpp:624
double * GetRowEntries(const int row)
Return a pointer to the entries in a row.
Definition: sparsemat.cpp:251
double * GetData() const
Definition: vector.hpp:88
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:903
FaceElementTransformations * GetFaceElementTransformations(int FaceNo, int mask=31)
Definition: mesh.cpp:557
void DeleteAll()
Delete whole array.
Definition: array.hpp:455
SparseMatrix * Ct
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
Definition: mesh.cpp:3755
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:185
void Finalize()
Finalize the construction of the hybridized matrix.
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
int Append(const T &el)
Append element to array, resize if necessary.
Definition: array.hpp:368
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:140
int GetSharedFace(int sface) const
Return the local face index for the given shared face.
Definition: pmesh.cpp:1422
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:3812
Array< int > hat_offsets
HypreParMatrix * RAP(HypreParMatrix *A, HypreParMatrix *P)
Returns the matrix P^t * A * P.
Definition: hypre.cpp:1318
void Reserve(int capacity)
Ensures that the allocated size is at least the given size.
Definition: array.hpp:122
void Transpose(const Table &A, Table &At, int _ncols_A)
Transpose a Table.
Definition: table.cpp:385
const SparseMatrix * GetConformingProlongation()
Definition: fespace.cpp:912
virtual void Finalize(int skip_zeros=1)
Definition: sparsemat.cpp:680
virtual const SparseMatrix * GetRestrictionMatrix()
Definition: fespace.hpp:150
void BooleanMult(int alpha, int *x, int beta, int *y)
Definition: hypre.hpp:363
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:1733
void Threshold(double eps)
Replace small entries, abs(a_ij) &lt;= eps, with zero.
Definition: densemat.cpp:2660
BilinearFormIntegrator * c_bfi
HypreParMatrix * Dof_TrueDof_Matrix()
The true dof-to-dof interpolation matrix.
Definition: pfespace.cpp:411
double MaxMaxNorm() const
Compute the norm ||A|| = max_{ij} |A_{ij}|.
Definition: densemat.cpp:691
double * Data() const
Returns vector of the elements.
Definition: densemat.hpp:77
int GetMyRank() const
Definition: pmesh.hpp:88
SparseMatrix * H
void AssembleMatrix(int el, const DenseMatrix &A)
Assemble the element matrix A into the hybridized system matrix.
Abstract finite element space.
Definition: fespace.hpp:62
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Definition: sparsemat.cpp:417
void GetIBDofs(int el, Array< int > &i_dofs, Array< int > &b_dofs) const
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:323
void PrintMatlab(std::ostream &out=std::cout) const
Prints matrix in matlab format.
Definition: sparsemat.cpp:2245
void SetSubVector(const Array< int > &dofs, const Vector &elemvect)
Definition: vector.cpp:493
~Hybridization()
Destructor.
void SetDataAndSize(double *d, int s)
Definition: vector.hpp:69
virtual void AssembleFaceMatrix(const FiniteElement &el1, const FiniteElement &el2, FaceElementTransformations &Trans, DenseMatrix &elmat)
Definition: bilininteg.cpp:39
FiniteElementSpace * fes
void MultTranspose(const Vector &x, Vector &y) const
Multiply a vector with the transposed matrix. y = At * x.
Definition: sparsemat.cpp:492
const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement associated with i&#39;th element.
Definition: fespace.cpp:1199
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:173
HypreParMatrix * pH
void MakeRef(T *, int)
Make this Array a reference to a pointer.
Definition: array.hpp:466
void MultAtB(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &AtB)
Multiply the transpose of a matrix A with a matrix B: At*B.
Definition: densemat.cpp:3362
Vector data type.
Definition: vector.hpp:33
void ReduceRHS(const Vector &b, Vector &b_r) const
void Init(const Array< int > &ess_tdof_list)
Prepare the Hybridization object for assembly.
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:143
Hybridization(FiniteElementSpace *fespace, FiniteElementSpace *c_fespace)
Constructor.
Class for parallel meshes.
Definition: pmesh.hpp:28
static void AdjustVDofs(Array< int > &vdofs)
Definition: fespace.cpp:148
double * data
Definition: densemat.hpp:375