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