MFEM  v3.3.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;
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_Int 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_Int)shared.slaves.size();
72  MPI_Allreduce(&num_shared_slave_faces, &glob_num_shared_slave_faces, 1,
73  HYPRE_MPI_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_Int Ct_num_rows = Ct->Height();
178  Array<HYPRE_Int> Ct_rows, *offsets[1] = { &Ct_rows };
179  pmesh->GenerateOffsets(1, &Ct_num_rows, offsets);
181  HYPRE_Int c_ldof_offset = c_pfes->GetMyDofOffset();
182  const HYPRE_Int *c_face_nbr_glob_ldof =
183  c_pfes->GetFaceNbrGlobalDofMap();
184  int *J = Ct->GetJ();
185  for (int i = 0; i < Ct_J.Size(); i++)
186  {
187  Ct_J[i] = J[i] < c_vsize ?
188  J[i] + c_ldof_offset :
189  c_face_nbr_glob_ldof[J[i] - c_vsize];
190  }
191  HypreParMatrix pCt(pmesh->GetComm(), Ct->Height(),
192  Ct_rows.Last(), c_pfes->GlobalVSize(),
193  Ct->GetI(), Ct_J.GetData(), Ct->GetData(),
194  Ct_rows, c_pfes->GetDofOffsets());
195  Ct_J.DeleteAll();
196  pC = pCt.Transpose();
197  }
198  if (pmesh->Nonconforming())
199  {
200  // TODO - Construct P_pc directly in the pH format
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  mfem::out << "\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  const bool fix_empty_rows = true;
658 #ifndef MFEM_USE_MPI
659  H->Finalize(skip_zeros, fix_empty_rows);
660 #else
661  ParFiniteElementSpace *c_pfes = dynamic_cast<ParFiniteElementSpace*>(c_fes);
662  if (!pC)
663  {
664  H->Finalize(skip_zeros, fix_empty_rows);
665  if (!c_pfes) { return; }
666 
667  OperatorHandle pP(pH.Type()), dH(pH.Type());
668  // TODO - construct P_pc / Dof_TrueDof_Matrix directly in the pH format
669  pP.ConvertFrom(P_pc ? P_pc : c_pfes->Dof_TrueDof_Matrix());
670  dH.MakeSquareBlockDiag(c_pfes->GetComm(),c_pfes->GlobalVSize(),
671  c_pfes->GetDofOffsets(), H);
672  pH.MakePtAP(dH, pP);
673  delete H;
674  H = NULL;
675  }
676  else
677  {
678  // TODO: add ones on the diagonal of zero rows
679  V->Finalize();
680  Array<HYPRE_Int> V_J(V->NumNonZeroElems());
681  MFEM_ASSERT(c_pfes, "");
682  const int c_vsize = c_fes->GetVSize();
683  HYPRE_Int c_ldof_offset = c_pfes->GetMyDofOffset();
684  const HYPRE_Int *c_face_nbr_glob_ldof = c_pfes->GetFaceNbrGlobalDofMap();
685  int *J = V->GetJ();
686  for (int i = 0; i < V_J.Size(); i++)
687  {
688  V_J[i] = J[i] < c_vsize ?
689  J[i] + c_ldof_offset :
690  c_face_nbr_glob_ldof[J[i] - c_vsize];
691  }
692  // TODO - lpH directly in the pH format
693  HypreParMatrix *lpH;
694  {
695  HypreParMatrix pV(c_pfes->GetComm(), V->Height(),
697  V->GetI(), V_J.GetData(), V->GetData(),
698  pC->ColPart(), pC->RowPart());
699  // The above constructor makes copies of all input arrays, so we can
700  // safely delete V_J and V:
701  V_J.DeleteAll();
702  delete V;
703  lpH = ParMult(pC, &pV);
704  }
705  OperatorHandle pP(pH.Type()), plpH(pH.Type());
706  // TODO - construct P_pc directly in the pH format
707  pP.ConvertFrom(P_pc);
708  plpH.ConvertFrom(lpH);
709  MFEM_VERIFY(pH.Type() != Operator::PETSC_MATIS, "To be implemented");
710  pH.MakePtAP(plpH, pP);
711  delete lpH;
712  }
713 #endif
714 }
715 
717 {
718 #ifndef MFEM_USE_MPI
719  if (!H) { ComputeH(); }
720 #else
721  if (!H && !pH.Ptr()) { ComputeH(); }
722 #endif
723 }
724 
725 void Hybridization::MultAfInv(const Vector &b, const Vector &lambda, Vector &bf,
726  int mode) const
727 {
728  // b1 = Rf^t b (assuming that Ref = 0)
729  Vector b1;
730  const SparseMatrix *R = fes->GetRestrictionMatrix();
731  if (!R)
732  {
733  b1.SetDataAndSize(b.GetData(), b.Size());
734  }
735  else
736  {
737  b1.SetSize(fes->GetVSize());
738  R->MultTranspose(b, b1);
739  }
740 
741  const int NE = fes->GetMesh()->GetNE();
742  Array<int> vdofs, i_dofs, b_dofs;
743  Vector el_vals, bf_i, i_vals, b_vals;
744  bf.SetSize(hat_offsets[NE]);
745  if (mode == 1)
746  {
747 #ifdef MFEM_USE_MPI
748  ParFiniteElementSpace *c_pfes =
749  dynamic_cast<ParFiniteElementSpace*>(c_fes);
750  if (!c_pfes)
751  {
752  Ct->Mult(lambda, bf);
753  }
754  else
755  {
756  Vector L(c_pfes->GetVSize());
757  (P_pc ? P_pc : c_pfes->GetProlongationMatrix())->Mult(lambda, L);
758  pC ? pC->MultTranspose(L, bf) : Ct->Mult(L, bf);
759  }
760 #else
761  Ct->Mult(lambda, bf);
762 #endif
763  }
764  // Apply Af^{-1}
765  Array<bool> vdof_marker(b1.Size());
766  vdof_marker = false;
767  for (int i = 0; i < NE; i++)
768  {
769  fes->GetElementVDofs(i, vdofs);
770  b1.GetSubVector(vdofs, el_vals);
771  for (int j = 0; j < vdofs.Size(); j++)
772  {
773  int vdof = vdofs[j];
774  if (vdof < 0) { vdof = -1 - vdof; }
775  if (vdof_marker[vdof]) { el_vals(j) = 0.0; }
776  else { vdof_marker[vdof] = true; }
777  }
778  bf_i.SetDataAndSize(&bf[hat_offsets[i]], vdofs.Size());
779  if (mode == 1)
780  {
781  el_vals -= bf_i;
782  }
783  GetIBDofs(i, i_dofs, b_dofs);
784  el_vals.GetSubVector(i_dofs, i_vals);
785  el_vals.GetSubVector(b_dofs, b_vals);
786 
787  LUFactors LU_ii(Af_data + Af_offsets[i], Af_ipiv + Af_f_offsets[i]);
788  double *U_ib = LU_ii.data + i_dofs.Size()*i_dofs.Size();
789  double *L_bi = U_ib + i_dofs.Size()*b_dofs.Size();
790  LUFactors LU_bb(L_bi + b_dofs.Size()*i_dofs.Size(),
791  LU_ii.ipiv + i_dofs.Size());
792  LU_ii.BlockForwSolve(i_dofs.Size(), b_dofs.Size(), 1, L_bi,
793  i_vals.GetData(), b_vals.GetData());
794  LU_bb.Solve(b_dofs.Size(), 1, b_vals.GetData());
795  bf_i = 0.0;
796  if (mode == 1)
797  {
798  LU_ii.BlockBackSolve(i_dofs.Size(), b_dofs.Size(), 1, U_ib,
799  b_vals.GetData(), i_vals.GetData());
800  bf_i.SetSubVector(i_dofs, i_vals);
801  }
802  bf_i.SetSubVector(b_dofs, b_vals);
803  }
804 }
805 
806 void Hybridization::ReduceRHS(const Vector &b, Vector &b_r) const
807 {
808  // bf = Af^{-1} Rf^t b
809  Vector bf;
810  MultAfInv(b, b, bf, 0);
811 
812  // b_r = Cf bf
813 #ifdef MFEM_USE_MPI
814  ParFiniteElementSpace *c_pfes = dynamic_cast<ParFiniteElementSpace*>(c_fes);
815  if (!c_pfes)
816  {
817  b_r.SetSize(Ct->Width());
818  Ct->MultTranspose(bf, b_r);
819  }
820  else
821  {
822  Vector bl(pC ? pC->Height() : Ct->Width());
823  pC ? pC->Mult(bf, bl) : Ct->MultTranspose(bf, bl);
824  b_r.SetSize(pH.Ptr()->Height());
825  (P_pc ? P_pc : c_pfes->GetProlongationMatrix())->MultTranspose(bl, b_r);
826  }
827 #else
828  b_r.SetSize(Ct->Width());
829  Ct->MultTranspose(bf, b_r);
830 #endif
831 }
832 
833 void Hybridization::ComputeSolution(const Vector &b, const Vector &sol_r,
834  Vector &sol) const
835 {
836  // bf = Af^{-1} ( Rf^t - Cf^t sol_r )
837  Vector bf;
838  MultAfInv(b, sol_r, bf, 1);
839 
840  // sol = Rf bf
841  GridFunction s;
842  const SparseMatrix *R = fes->GetRestrictionMatrix();
843  if (!R)
844  {
845  MFEM_ASSERT(sol.Size() == fes->GetVSize(), "");
846  s.MakeRef(fes, sol, 0);
847  }
848  else
849  {
850  s.SetSpace(fes);
851  R->MultTranspose(sol, s);
852  }
853  const int NE = fes->GetMesh()->GetNE();
854  Array<int> vdofs;
855  for (int i = 0; i < NE; i++)
856  {
857  fes->GetElementVDofs(i, vdofs);
858  for (int j = hat_offsets[i]; j < hat_offsets[i+1]; j++)
859  {
860  if (hat_dofs_marker[j] == 1) { continue; } // skip essential b.c.
861  int vdof = vdofs[j-hat_offsets[i]];
862  if (vdof >= 0) { s(vdof) = bf(j); }
863  else { s(-1-vdof) = -bf(j); }
864  }
865  }
866  if (R)
867  {
868  R->Mult(s, sol); // assuming that Ref = 0
869  }
870 }
871 
873 {
874  delete H;
875  H = NULL;
876 #ifdef MFEM_USE_MPI
877  pH.Clear();
878 #endif
879 }
880 
881 }
Abstract class for Finite Elements.
Definition: fe.hpp:47
virtual const Operator * GetProlongationMatrix()
Definition: pfespace.cpp:632
int Size() const
Logical size of the array.
Definition: array.hpp:110
Array< int > Af_offsets
int RowSize(const int i) const
Returns the number of elements in row i.
Definition: sparsemat.cpp:213
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:498
int GetVSize() const
Definition: fespace.hpp:163
virtual int NumNonZeroElems() const
Returns the number of the nonzero elements in the matrix.
Definition: sparsemat.cpp:972
void ConvertFrom(OperatorHandle &A)
Convert the given OperatorHandle A to the currently set type id.
Definition: handle.cpp:174
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:27
const SparseMatrix * GetConformingProlongation() const
Definition: fespace.cpp:725
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:353
HYPRE_Int * GetDofOffsets() const
Definition: pfespace.hpp:174
virtual void Finalize(int skip_zeros=1)
Finalize the matrix initialization, switching the storage format from LIL to CSR. ...
Definition: sparsemat.hpp:278
HYPRE_Int * ColPart()
Returns the column partitioning.
Definition: hypre.hpp:334
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:1017
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:320
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:133
HypreParMatrix * P_pc
int * GetRowColumns(const int row)
Return a pointer to the column indices in a row.
Definition: sparsemat.cpp:260
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:633
Lists all edges/faces in the nonconforming mesh.
Definition: ncmesh.hpp:166
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:468
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:42
int GetNSharedFaces() const
Return the number of shared faces (3D), edges (2D), vertices (1D)
Definition: pmesh.cpp:2009
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
Operator::Type Type() const
Get the currently set operator type id.
Definition: handle.hpp:85
void BlockFactor(int m, int n, double *A12, double *A21, double *A22) const
Definition: densemat.cpp:4077
void BlockBackSolve(int m, int n, int r, const double *U12, const double *X2, double *Y1) const
Definition: densemat.cpp:4113
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Definition: vector.cpp:462
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
int Size() const
Returns the size of the vector.
Definition: vector.hpp:113
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:614
Abstract parallel finite element space.
Definition: pfespace.hpp:31
const FiniteElement * GetFaceElement(int i) const
Definition: fespace.cpp:1432
void Factor(int m)
Definition: densemat.cpp:3832
HypreParMatrix * GetPartialConformingInterpolation()
For a non-conforming mesh, construct and return the interpolation matrix from the partially conformin...
Definition: pfespace.cpp:1570
FaceElementTransformations * GetInteriorFaceTransformations(int FaceNo)
Definition: mesh.hpp:822
double * GetRowEntries(const int row)
Return a pointer to the entries in a row.
Definition: sparsemat.cpp:274
double * GetData() const
Definition: vector.hpp:121
const HYPRE_Int * GetFaceNbrGlobalDofMap()
Definition: pfespace.hpp:260
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:947
FaceElementTransformations * GetFaceElementTransformations(int FaceNo, int mask=31)
Definition: mesh.cpp:633
HypreParMatrix * ParMult(const HypreParMatrix *A, const HypreParMatrix *B)
Returns the matrix A * B.
Definition: hypre.cpp:1456
int GetConformingVSize() const
Definition: fespace.hpp:172
HYPRE_Int GetMyDofOffset() const
Definition: pfespace.cpp:622
void Reset()
Destroy the current hybridization matrix while preserving the computed constraint matrix and the set ...
void DeleteAll()
Delete whole array.
Definition: array.hpp:633
SparseMatrix * Ct
ParNCMesh * pncmesh
Definition: pmesh.hpp:134
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
Definition: mesh.cpp:3334
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:184
void Finalize()
Finalize the construction of the hybridized matrix.
bool Nonconforming() const
Definition: mesh.hpp:999
HYPRE_Int GetGlobalNumRows() const
Definition: hypre.hpp:374
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.
MPI_Comm GetComm() const
Definition: pfespace.hpp:166
HYPRE_Int GetGlobalNumCols() const
Definition: hypre.hpp:377
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:36
int Append(const T &el)
Append element to array, resize if necessary.
Definition: array.hpp:548
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:2028
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:4104
Array< int > hat_offsets
FaceElementTransformations * GetSharedFaceTransformations(int sf, bool fill2=true)
Definition: pmesh.cpp:1924
HypreParMatrix * Dof_TrueDof_Matrix() const
The true dof-to-dof interpolation matrix.
Definition: pfespace.hpp:198
void Reserve(int capacity)
Ensures that the allocated size is at least the given size.
Definition: array.hpp:123
void Transpose(const Table &A, Table &At, int _ncols_A)
Transpose a Table.
Definition: table.cpp:408
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:3969
HypreParMatrix * Transpose() const
Returns the transpose of *this.
Definition: hypre.cpp:936
virtual void MakeRef(FiniteElementSpace *f, double *v)
Make the GridFunction reference external data on a new FiniteElementSpace.
Definition: gridfunc.cpp:181
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:153
void PrintMatlab(std::ostream &out=mfem::out) const
Prints matrix in matlab format.
Definition: sparsemat.cpp:2382
virtual const SparseMatrix * GetRestrictionMatrix()
Definition: fespace.hpp:149
int Dimension() const
Definition: mesh.hpp:641
void BooleanMult(int alpha, int *x, int beta, int *y)
Definition: hypre.hpp:404
void MultAfInv(const Vector &b, const Vector &lambda, Vector &bf, int mode) const
double * GetData() const
Return element data, i.e. array A.
Definition: sparsemat.hpp:135
void AddSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Definition: sparsemat.cpp:1850
Operator * Ptr() const
Access the underlying Operator pointer.
Definition: handle.hpp:82
int * GetI() const
Return the array I.
Definition: sparsemat.hpp:131
void Threshold(double eps)
Replace small entries, abs(a_ij) &lt;= eps, with zero.
Definition: densemat.cpp:2833
HYPRE_Int * RowPart()
Returns the row partitioning.
Definition: hypre.hpp:332
BilinearFormIntegrator * c_bfi
HYPRE_Int GlobalVSize() const
Definition: pfespace.hpp:176
double MaxMaxNorm() const
Compute the norm ||A|| = max_{ij} |A_{ij}|.
Definition: densemat.cpp:795
double * Data() const
Returns the matrix data array.
Definition: densemat.hpp:92
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:117
std::vector< Slave > slaves
Definition: ncmesh.hpp:170
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Definition: sparsemat.cpp:475
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:503
OperatorHandle pH
void Clear()
Clear the OperatorHandle, deleting the held Operator (if owned), while leaving the type id unchanged...
Definition: handle.hpp:110
~Hybridization()
Destructor.
void SetDataAndSize(double *d, int s)
Set the Vector data and size.
Definition: vector.hpp:94
int GetElementBaseGeometry(int i=0) const
Definition: mesh.hpp:684
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:552
const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement associated with i&#39;th element.
Definition: fespace.cpp:1171
const FiniteElement * GetFaceNbrFaceFE(int i) const
Definition: pfespace.cpp:913
virtual void SetSpace(FiniteElementSpace *f)
Associate a new FiniteElementSpace with the GridFunction.
Definition: gridfunc.cpp:173
int GetFaceNbrVSize() const
Definition: pfespace.hpp:255
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:145
T & Last()
Return the last element in the array.
Definition: array.hpp:581
void GetFaceNbrFaceVDofs(int i, Array< int > &vdofs) const
Definition: pfespace.cpp:871
void MakeRef(T *, int)
Make this Array a reference to a pointer.
Definition: array.hpp:644
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:3637
Vector data type.
Definition: vector.hpp:41
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:177
void GenerateOffsets(int N, HYPRE_Int loc_sizes[], Array< HYPRE_Int > *offsets[]) const
Definition: pmesh.cpp:1241
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:139
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:86
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:64
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:175
int * GetJ() const
Return the array J.
Definition: sparsemat.hpp:133
Hybridization(FiniteElementSpace *fespace, FiniteElementSpace *c_fespace)
Constructor.
Class for parallel meshes.
Definition: pmesh.hpp:29
void AdjustDofDirection(Array< int > &dofs)
Definition: densemat.cpp:2774
void MakePtAP(OperatorHandle &A, OperatorHandle &P)
Reset the OperatorHandle to hold the product P^t A P.
Definition: handle.cpp:109
static void AdjustVDofs(Array< int > &vdofs)
Definition: fespace.cpp:120
void SetType(Operator::Type tid)
Invoke Clear() and set a new type id.
Definition: handle.hpp:118
double * data
Definition: densemat.hpp:437