MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
pweakform.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11
12#include "pweakform.hpp"
13
14#ifdef MFEM_USE_MPI
15
16namespace mfem
17{
18
20{
21 int j;
22 for (int i = 0; i<ess_tdof_list.Size(); i++)
23 {
24 int tdof = ess_tdof_list[i];
25 for (j = 0; j < nblocks; j++)
26 {
27 if (tdof_offsets[j+1] > tdof) { break; }
28 }
29 ess_tdofs[j]->Append(tdof-tdof_offsets[j]);
30 }
31}
32
33void ParDPGWeakForm::Assemble(int skip_zeros)
34{
35 DPGWeakForm::Assemble(skip_zeros);
36}
37
39{
40 if (!P) { BuildProlongation(); }
41
44 p_mat->owns_blocks = 1;
46 HypreParMatrix * A = nullptr;
47 HypreParMatrix * PtAP = nullptr;
48 for (int i = 0; i<nblocks; i++)
49 {
50 HypreParMatrix * Pi = (HypreParMatrix*)(&P->GetBlock(i,i));
51 HypreParMatrix * Pit = Pi->Transpose();
52 for (int j = 0; j<nblocks; j++)
53 {
54 if (m->IsZeroBlock(i,j)) { continue; }
55 if (i == j)
56 {
57 // Make block diagonal square hypre matrix
58 A = new HypreParMatrix(trial_pfes[i]->GetComm(), trial_pfes[i]->GlobalVSize(),
59 trial_pfes[i]->GetDofOffsets(),&m->GetBlock(i,i));
60 PtAP = RAP(A,Pi);
61 delete A;
63 }
64 else
65 {
66 HypreParMatrix * Pj = (HypreParMatrix*)(&P->GetBlock(j,j));
67 A = new HypreParMatrix(trial_pfes[i]->GetComm(), trial_pfes[i]->GlobalVSize(),
68 trial_pfes[j]->GlobalVSize(), trial_pfes[i]->GetDofOffsets(),
69 trial_pfes[j]->GetDofOffsets(), &m->GetBlock(i,j));
70 HypreParMatrix * APj = ParMult(A, Pj,true);
71 delete A;
72 PtAP = ParMult(Pit,APj,true);
73 delete APj;
74 p_mat_e->SetBlock(i,j,PtAP->EliminateCols(*ess_tdofs[j]));
75 PtAP->EliminateRows(*ess_tdofs[i]);
76 }
77 p_mat->SetBlock(i,j,PtAP);
78 }
79 delete Pit;
80 }
81}
82
83
85{
88 P->owns_blocks = 0;
89 R->owns_blocks = 0;
90
91 for (int i = 0; i<nblocks; i++)
92 {
93 HypreParMatrix * P_ = trial_pfes[i]->Dof_TrueDof_Matrix();
94 P->SetBlock(i,i,P_);
95 const SparseMatrix * R_ = trial_pfes[i]->GetRestrictionMatrix();
96 R->SetBlock(i,i,const_cast<SparseMatrix*>(R_));
97 }
98}
99
101 &ess_tdof_list,
102 Vector &x,
103 OperatorHandle &A, Vector &X,
104 Vector &B, int copy_interior)
105{
106 FormSystemMatrix(ess_tdof_list, A);
107
108
109 if (static_cond)
110 {
111 static_cond->ReduceSystem(x, X, B, copy_interior);
112 }
113 else
114 {
115 B.SetSize(P->Width());
116 P->MultTranspose(*y,B);
117 X.SetSize(R->Height());
118 R->Mult(x,X);
119
120 // eliminate tdof in RHS
121 // B -= Ae*X
122 Vector tmp(B.Size());
123 p_mat_e->Mult(X,tmp);
124 B-=tmp;
125
126 for (int j = 0; j<nblocks; j++)
127 {
128 if (!ess_tdofs[j]->Size()) { continue; }
129 for (int i = 0; i < ess_tdofs[j]->Size(); i++)
130 {
131 int tdof = (*ess_tdofs[j])[i];
132 int gdof = tdof + tdof_offsets[j];
133 B(gdof) = X(gdof); // diagonal policy in always one in parallel
134 }
135 }
136 if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
137 }
138}
139
141 &ess_tdof_list,
143{
144 if (static_cond)
145 {
147 {
148 static_cond->SetEssentialTrueDofs(ess_tdof_list);
150 }
152 }
153 else
154 {
155 FillEssTdofLists(ess_tdof_list);
156 if (mat)
157 {
158 const int remove_zeros = 0;
159 Finalize(remove_zeros);
161 delete mat;
162 mat = nullptr;
163 delete mat_e;
164 mat_e = nullptr;
165 }
166 A.Reset(p_mat,false);
167 }
168
169}
170
171
172
174 Vector &x)
175{
176
177 if (static_cond)
178 {
180 }
181 else
182 {
183 x.SetSize(P->Height());
184 P->Mult(X, x);
185 }
186}
187
189{
191 delete p_mat_e;
192 p_mat_e = nullptr;
193 delete p_mat;
194 p_mat = nullptr;
195 for (int i = 0; i<nblocks; i++)
196 {
197 delete ess_tdofs[i];
198 ess_tdofs[i] = new Array<int>();
199 }
200 delete P;
201 P = nullptr;
202 delete R;
203 R = nullptr;
204}
205
207{
208 delete p_mat_e;
209 p_mat_e = nullptr;
210 delete p_mat;
211 p_mat = nullptr;
212 for (int i = 0; i<nblocks; i++)
213 {
214 delete ess_tdofs[i];
215 }
216 delete P;
217 delete R;
218}
219
220} // namespace mfem
221
222#endif
int Size() const
Return the logical size of the array.
Definition array.hpp:144
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Definition array.hpp:769
void SetBlock(int i, int j, SparseMatrix *mat)
Set A(i,j) = mat.
int IsZeroBlock(int i, int j) const
Check if block (i,j) is a zero block.
SparseMatrix & GetBlock(int i, int j)
Return a reference to block (i,j). Reference may be invalid if Aij(i,j) == NULL.
int owns_blocks
If owns_blocks the SparseMatrix objects Aij will be deallocated.
virtual void Mult(const Vector &x, Vector &y) const
Matrix-Vector Multiplication y = A*x.
A class to handle Block systems in a matrix-free implementation.
void SetBlock(int iRow, int iCol, Operator *op, real_t c=1.0)
Add a block op in the block-entry (iblock, jblock).
Operator & GetBlock(int i, int j)
Return a reference to block i,j.
virtual void Mult(const Vector &x, Vector &y) const
Operator application.
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator.
void FormSystemMatrix(Operator::DiagonalPolicy diag_policy)
void SetEssentialTrueDofs(const Array< int > &ess_tdof_list)
Determine and save internally essential reduced true dofs.
BlockOperator & GetParallelSchurMatrix()
Return the parallel Schur complement matrix.
void ReduceSystem(Vector &x, Vector &X, Vector &B, int copy_interior=0) const
Set the reduced solution X and r.h.s B vectors from the full linear system solution x and r....
void ComputeSolution(const Vector &sc_sol, Vector &sol) const
BlockStaticCondensation * static_cond
Owned.
Definition weakform.hpp:38
BlockMatrix * mat
Block matrix to be associated with the Block bilinear form. Owned.
Definition weakform.hpp:49
Array< int > tdof_offsets
Definition weakform.hpp:46
virtual void Update()
Update the DPGWeakForm after mesh modifications (AMR)
Definition weakform.cpp:671
Array< int > dof_offsets
Definition weakform.hpp:45
int Size() const
Get the size of the bilinear form of the DPGWeakForm.
Definition weakform.hpp:147
BlockVector * y
Block vector to be associated with the Block linear form.
Definition weakform.hpp:52
void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
Definition weakform.cpp:97
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all integrators.
Definition weakform.cpp:247
BlockMatrix * mat_e
Block Matrix used to store the eliminations from the b.c. Owned. .
Definition weakform.hpp:57
Wrapper for hypre's ParCSR matrix class.
Definition hypre.hpp:388
void EliminateRows(const Array< int > &rows)
Eliminate rows from the diagonal and off-diagonal blocks of the matrix.
Definition hypre.cpp:2395
void EliminateRowsCols(const Array< int > &rows_cols, const HypreParVector &X, HypreParVector &B)
Definition hypre.cpp:2356
HypreParMatrix * Transpose() const
Returns the transpose of *this.
Definition hypre.cpp:1684
HypreParMatrix * EliminateCols(const Array< int > &cols)
Definition hypre.cpp:2381
Pointer to an Operator of a specified type.
Definition handle.hpp:34
void Reset(OpType *A, bool own_A=true)
Reset the OperatorHandle to the given OpType pointer, A.
Definition handle.hpp:145
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition operator.hpp:66
@ DIAG_ONE
Set the diagonal value to one.
Definition operator.hpp:50
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition operator.hpp:72
BlockOperator * p_mat
Definition pweakform.hpp:44
Array< ParFiniteElementSpace * > trial_pfes
Definition pweakform.hpp:30
BlockOperator * P
Definition pweakform.hpp:40
virtual ~ParDPGWeakForm()
Destroys bilinear form.
void Assemble(int skip_zeros=1)
Assemble the local matrix.
Definition pweakform.cpp:33
void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
Form the linear system A X = B, corresponding to this DPG weak form.
BlockMatrix * R
Definition pweakform.hpp:41
virtual void RecoverFEMSolution(const Vector &X, Vector &x)
void FillEssTdofLists(const Array< int > &ess_tdof_list)
Definition pweakform.cpp:19
Array< Array< int > * > ess_tdofs
Definition pweakform.hpp:33
BlockOperator * p_mat_e
Definition pweakform.hpp:45
void ParallelAssemble(BlockMatrix *mat)
Returns the matrix assembled on the true dofs, i.e. P^t A P.
Definition pweakform.cpp:38
void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Form the linear system matrix A, see FormLinearSystem() for details.
virtual void Update()
Update the DPGWeakForm after mesh modifications (AMR)
Data type sparse matrix.
Definition sparsemat.hpp:51
Vector data type.
Definition vector.hpp:80
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
void SetSubVectorComplement(const Array< int > &dofs, const real_t val)
Set all vector entries NOT in the dofs Array to the given val.
Definition vector.cpp:739
void RAP(const DenseMatrix &A, const DenseMatrix &P, DenseMatrix &RAP)
HypreParMatrix * ParMult(const HypreParMatrix *A, const HypreParMatrix *B, bool own_matrix)
Definition hypre.cpp:2949