23 #ifdef MFEM_DEBUG_HYBRIDIZATION_CP
32 : fes(fespace), c_fes(c_fespace), c_bfi(NULL), Ct(NULL), H(NULL),
33 Af_data(NULL), Af_ipiv(NULL)
61 int c_num_face_nbr_dofs = 0;
65 HYPRE_Int num_shared_slave_faces = 0, glob_num_shared_slave_faces = 0;
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)
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
89 #undef MFEM_DEBUG_HERE
96 Ct =
new SparseMatrix(num_hat_dofs, c_vsize + c_num_face_nbr_dofs);
100 const int skip_zeros = 1;
105 for (
int i = 0; i < num_faces; i++)
108 if (!FTr) {
continue; }
115 for (
int j = 0; j < s1; j++)
119 for (
int j = 0; j < s2; j++)
121 vdofs[s1+j] = o2 + j;
137 for (
int i = 0; i < num_shared_faces; i++)
140 const bool ghost_sface = (face_no >= num_faces);
145 MFEM_ASSERT(FTr->
Elem2No < 0,
"");
151 const int fill2 =
false;
156 for (
int j = 0; j < c_vdofs.
Size(); j++)
158 c_vdofs[j] += c_vsize;
164 for (
int j = 0; j < s1; j++)
174 if (glob_num_shared_slave_faces)
178 HYPRE_Int Ct_num_rows =
Ct->
Height();
183 const HYPRE_Int *c_face_nbr_glob_ldof =
186 for (
int i = 0; i < Ct_J.Size(); i++)
188 Ct_J[i] = J[i] < c_vsize ?
189 J[i] + c_ldof_offset :
190 c_face_nbr_glob_ldof[J[i] - c_vsize];
210 MFEM_ABORT(
"TODO: algebraic definition of C");
221 int num_hat_dofs = 0;
224 for (
int i = 0; i < NE; i++)
227 num_hat_dofs += vdofs.
Size();
234 #ifdef MFEM_DEBUG_HYBRIDIZATION_CP
237 std::ofstream C_file(
"C_matrix.txt");
245 std::ofstream P_file(
"P_matrix.txt");
263 free_tdof_marker = 1;
264 for (
int i = 0; i < ess_tdof_list.
Size(); i++)
266 free_tdof_marker[ess_tdof_list[i]] = 0;
275 free_vdofs_marker.
MakeRef(free_tdof_marker);
280 cP->
BooleanMult(free_tdof_marker, free_vdofs_marker);
287 P->
BooleanMult(1, free_tdof_marker, 0, free_vdofs_marker);
293 free_vdofs_marker.
MakeRef(free_tdof_marker);
298 cP->
BooleanMult(free_tdof_marker, free_vdofs_marker);
301 for (
int i = 0; i < NE; i++)
305 for (
int j = 0; j < vdofs.
Size(); j++)
319 for (
int i = 0; i < num_hat_dofs; i++)
334 #ifdef MFEM_DEBUG_HERE
337 for (
int i = 0; i < NE; i++)
343 #ifdef MFEM_DEBUG_HERE
351 #ifdef MFEM_DEBUG_HERE
355 int myid = pmesh ? pmesh->GetMyRank() : 0;
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
375 for (
int i = 0; i < NE; i++)
379 for (
int j = 0; j < vdofs.
Size(); j++)
383 vdof_marker[vdofs[j]] = 1;
387 for (
int tdof = 0; tdof < R->
Height(); tdof++)
389 if (free_tdof_marker[tdof]) {
continue; }
391 const int ncols = R->
RowSize(tdof);
394 for (
int j = 0; j < ncols; j++)
396 if (std::abs(vals[j]) != 0.0 && vdof_marker[cols[j]] == 0)
398 MFEM_ABORT(
"Ref != 0");
417 for (
int i = h_start; i < h_end; i++)
420 if (mark == 0) { i_dofs.
Append(i-h_start); }
421 else if (mark == -1) { b_dofs.
Append(i-h_start); }
433 for (
int i = h_start; i < h_end; i++)
436 if (mark == 0) { num_idofs++; }
437 else if (mark == -1) { b_dofs.
Append(i); }
455 for (
int j = 0; j < i_dofs.
Size(); j++)
457 int j_dof = i_dofs[j];
458 for (
int i = 0; i < i_dofs.
Size(); i++)
460 A_ii(i,j) = A(i_dofs[i],j_dof);
462 for (
int i = 0; i < b_dofs.
Size(); i++)
464 A_bi(i,j) = A(b_dofs[i],j_dof);
467 for (
int j = 0; j < b_dofs.
Size(); j++)
469 int j_dof = b_dofs[j];
470 for (
int i = 0; i < i_dofs.
Size(); i++)
472 A_ib(i,j) = A(i_dofs[i],j_dof);
474 for (
int i = 0; i < b_dofs.
Size(); i++)
476 A_bb(i,j) = A(b_dofs[i],j_dof);
503 Ordering::DofsToVDofs<Ordering::byNODES>(e2f.
Size()/vdim, vdim, lvdofs);
504 MFEM_ASSERT(lvdofs.
Size() == A.
Height(),
"internal error");
507 for (
int i = 0; i < lvdofs.
Size(); i++)
510 bd = (bd >= 0) ? vdofs[bd] : -1-vdofs[-1-bd];
511 MFEM_ASSERT(bvdofs[i] == bd,
"internal error");
517 for (
int i = 0; i < lvdofs.
Size(); i++)
533 for (
int j = 0; j < i_dofs.
Size(); j++)
535 int j_f = e2f[i_dofs[j]];
536 if (j_f == -1) {
continue; }
537 for (
int i = 0; i < i_dofs.Size(); i++)
539 int i_f = e2f[i_dofs[i]];
540 if (i_f == -1) {
continue; }
541 A_ii(i,j) += B(i_f,j_f);
543 for (
int i = 0; i < b_dofs.
Size(); i++)
545 int i_f = e2f[b_dofs[i]];
546 if (i_f == -1) {
continue; }
547 A_bi(i,j) += B(i_f,j_f);
550 for (
int j = 0; j < b_dofs.
Size(); j++)
552 int j_f = e2f[b_dofs[j]];
553 if (j_f == -1) {
continue; }
554 for (
int i = 0; i < i_dofs.
Size(); i++)
556 int i_f = e2f[i_dofs[i]];
557 if (i_f == -1) {
continue; }
558 A_ib(i,j) += B(i_f,j_f);
560 for (
int i = 0; i < b_dofs.Size(); i++)
562 int i_f = e2f[b_dofs[i]];
563 if (i_f == -1) {
continue; }
564 A_bb(i,j) += B(i_f,j_f);
571 const int skip_zeros = 1;
585 int c_mark_start = 0;
586 for (
int el = 0; el < NE; 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();
595 LU_ii.
ipiv + i_dofs_size);
597 LU_ii.
Factor(i_dofs_size);
599 A_ib_data, A_bi_data, LU_bb.data);
600 LU_bb.Factor(b_dofs.
Size());
604 for (
int i = 0; i < b_dofs.
Size(); i++)
606 const int row = b_dofs[i];
609 for (
int j = 0; j < ncols; j++)
611 const int c_dof = cols[j];
612 if (c_dof_marker[c_dof] < c_mark_start)
614 c_dof_marker[c_dof] = c_mark_start + c_dofs.
Size();
621 for (
int i = 0; i < b_dofs.
Size(); i++)
623 const int row = b_dofs[i];
627 for (
int j = 0; j < ncols; j++)
629 const int loc_j = c_dof_marker[cols[j]] - c_mark_start;
630 Cb_t(i,loc_j) = vals[j];
642 MultAtB(Cb_t, Sb_inv_Cb_t, Hb);
650 V->AddSubMatrix(b_dofs, c_dofs, Sb_inv_Cb_t, skip_zeros);
654 c_mark_start += c_dofs.
Size();
655 MFEM_VERIFY(c_mark_start >= 0,
"overflow");
664 if (!c_pfes) {
return; }
678 MFEM_ASSERT(c_pfes,
"");
683 for (
int i = 0; i < V_J.Size(); i++)
685 V_J[i] = J[i] < c_vsize ?
686 J[i] + c_ldof_offset :
687 c_face_nbr_glob_ldof[J[i] - c_vsize];
693 V->GetI(), V_J.GetData(), V->GetData(),
733 Vector el_vals, bf_i, i_vals, b_vals;
757 for (
int i = 0; i < NE; i++)
761 for (
int j = 0; j < vdofs.
Size(); j++)
764 if (vdof < 0) { vdof = -1 - vdof; }
765 if (vdof_marker[vdof]) { el_vals(j) = 0.0; }
766 else { vdof_marker[vdof] =
true; }
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();
845 for (
int i = 0; i < NE; i++)
851 int vdof = vdofs[j-hat_offsets[i]];
852 if (vdof >= 0) { s(vdof) = bf(j); }
853 else { s(-1-vdof) = -bf(j); }
Abstract class for Finite Elements.
int Size() const
Logical size of the array.
int RowSize(const int i) const
Returns the number of elements in row i.
virtual int NumNonZeroElems() const
Returns the number of the nonzero elements in the matrix.
Class for grid function - Vector with associated FE space.
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.
void MakeRef(FiniteElementSpace *f, Vector &v, int v_offset)
HYPRE_Int * ColPart()
Returns the column partitioning.
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.
void SetSize(int s)
Resize the vector if the new size is different.
void GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i'th element.
int * GetRowColumns(const int row)
Return a pointer to the column indices in a row.
void BooleanMult(const Array< int > &x, Array< int > &y) const
y = A * x, but treat all elements as booleans (zero=false, nonzero=true).
Lists all edges/faces in the nonconforming mesh.
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols.
int GetNSharedFaces() const
Return the number of shared faces (3D), edges (2D), vertices (1D)
void BlockFactor(int m, int n, double *A12, double *A21, double *A22) const
void BlockBackSolve(int m, int n, int r, const double *U12, const double *X2, double *Y1) const
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
HYPRE_Int * GetDofOffsets()
Data type dense matrix using column-major storage.
int Size() const
Returns the size of the vector.
int GetNE() const
Returns number of elements.
Abstract parallel finite element space.
const FiniteElement * GetFaceElement(int i) const
HypreParMatrix * GetPartialConformingInterpolation()
For a non-conforming mesh, construct and return the interpolation matrix from the partially conformin...
FaceElementTransformations * GetInteriorFaceTransformations(int FaceNo)
double * GetRowEntries(const int row)
Return a pointer to the entries in a row.
const HYPRE_Int * GetFaceNbrGlobalDofMap()
HYPRE_Int Mult(HypreParVector &x, HypreParVector &y, double alpha=1.0, double beta=0.0)
Computes y = alpha * A * x + beta * y.
FaceElementTransformations * GetFaceElementTransformations(int FaceNo, int mask=31)
HYPRE_Int GetMyDofOffset() const
void Reset()
Destroy the current hybridization matrix while preserving the computed constraint matrix and the set ...
void DeleteAll()
Delete whole array.
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
int GetNE() const
Returns number of elements in the mesh.
void Finalize()
Finalize the construction of the hybridized matrix.
bool Nonconforming() const
HYPRE_Int GetGlobalNumRows() const
void ExchangeFaceNbrData()
void AssembleBdrMatrix(int bdr_el, const DenseMatrix &A)
Assemble the boundary element matrix A into the hybridized system matrix.
HYPRE_Int GetGlobalNumCols() const
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows.
int Append(const T &el)
Append element to array, resize if necessary.
Mesh * GetMesh() const
Returns the mesh.
int GetSharedFace(int sface) const
Return the local face index for the given shared face.
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
FaceElementTransformations * GetSharedFaceTransformations(int sf, bool fill2=true)
HypreParMatrix * RAP(HypreParMatrix *A, HypreParMatrix *P)
Returns the matrix P^t * A * P.
void Reserve(int capacity)
Ensures that the allocated size is at least the given size.
void Transpose(const Table &A, Table &At, int _ncols_A)
Transpose a Table.
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...
const SparseMatrix * GetConformingProlongation()
int GetVDim() const
Returns vector dimension.
virtual void Finalize(int skip_zeros=1)
virtual const SparseMatrix * GetRestrictionMatrix()
void BooleanMult(int alpha, int *x, int beta, int *y)
void MultAfInv(const Vector &b, const Vector &lambda, Vector &bf, int mode) const
double * GetData() const
Return element data.
void AddSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
int * GetI() const
Return the array I.
void Threshold(double eps)
Replace small entries, abs(a_ij) <= eps, with zero.
HYPRE_Int * RowPart()
Returns the row partitioning.
BilinearFormIntegrator * c_bfi
HypreParMatrix * Dof_TrueDof_Matrix()
The true dof-to-dof interpolation matrix.
double MaxMaxNorm() const
Compute the norm ||A|| = max_{ij} |A_{ij}|.
double * Data() const
Returns vector of the elements.
void AssembleMatrix(int el, const DenseMatrix &A)
Assemble the element matrix A into the hybridized system matrix.
std::vector< Slave > slaves
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
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.
void PrintMatlab(std::ostream &out=std::cout) const
Prints matrix in matlab format.
void SetSubVector(const Array< int > &dofs, const Vector &elemvect)
~Hybridization()
Destructor.
void SetDataAndSize(double *d, int s)
HypreParMatrix * Transpose()
Returns the transpose of *this.
int GetElementBaseGeometry(int i=0) const
void MultTranspose(const Vector &x, Vector &y) const
Multiply a vector with the transposed matrix. y = At * x.
const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement associated with i'th element.
const FiniteElement * GetFaceNbrFaceFE(int i) const
void SetSpace(FiniteElementSpace *f)
int GetFaceNbrVSize() const
Array< int > hat_dofs_marker
void GetFaceVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i'th face element (2D and 3D).
T & Last()
Return the last element in the array.
void GetFaceNbrFaceVDofs(int i, Array< int > &vdofs) const
void MakeRef(T *, int)
Make this Array a reference to a pointer.
const NCList & GetSharedList(int type)
Helper to get shared vertices/edges/faces ('type' == 0/1/2 resp.).
void MultAtB(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &AtB)
Multiply the transpose of a matrix A with a matrix B: At*B.
const FiniteElementCollection * FEColl() const
void GenerateOffsets(int N, HYPRE_Int loc_sizes[], Array< HYPRE_Int > *offsets[]) const
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'th boundary element.
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Wrapper for hypre's ParCSR matrix class.
int * GetJ() const
Return the array J.
Hybridization(FiniteElementSpace *fespace, FiniteElementSpace *c_fespace)
Constructor.
Class for parallel meshes.
void AdjustDofDirection(Array< int > &dofs)
static void AdjustVDofs(Array< int > &vdofs)
HypreParMatrix * ParMult(HypreParMatrix *A, HypreParMatrix *B)
Returns the matrix A * B.