MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
pcomplexweakform.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 "pcomplexweakform.hpp"
13
14#ifdef MFEM_USE_MPI
15
16namespace mfem
17{
18
20 ess_tdof_list)
21{
22 int j;
23 for (int i = 0; i < ess_tdof_list.Size(); i++)
24 {
25 int tdof = ess_tdof_list[i];
26 for (j = 0; j < nblocks; j++)
27 {
28 if (tdof_offsets[j+1] > tdof) { break; }
29 }
30 ess_tdofs[j]->Append(tdof-tdof_offsets[j]);
31 }
32}
33
35{
37}
38
40 BlockMatrix *m_i)
41{
42 if (!P) { BuildProlongation(); }
43
52 HypreParMatrix * A_r = nullptr;
53 HypreParMatrix * A_i = nullptr;
54 HypreParMatrix * PtAP_r = nullptr;
55 HypreParMatrix * PtAP_i = nullptr;
56 for (int i = 0; i < nblocks; i++)
57 {
58 HypreParMatrix * Pi = (HypreParMatrix*)(&P->GetBlock(i,i));
59 for (int j = 0; j<nblocks; j++)
60 {
61 if (m_r->IsZeroBlock(i,j)) { continue; }
62 if (i == j)
63 {
64 // Make block diagonal square hypre matrix
65 A_r = new HypreParMatrix(trial_pfes[i]->GetComm(), trial_pfes[i]->GlobalVSize(),
66 trial_pfes[i]->GetDofOffsets(), &m_r->GetBlock(i,i));
67 PtAP_r = RAP(A_r,Pi);
68 delete A_r;
69 p_mat_e_r->SetBlock(i, i, PtAP_r->EliminateRowsCols(*ess_tdofs[i]));
70
71 A_i = new HypreParMatrix(trial_pfes[i]->GetComm(), trial_pfes[i]->GlobalVSize(),
72 trial_pfes[i]->GetDofOffsets(), &m_i->GetBlock(i,i));
73
74 PtAP_i = RAP(A_i,Pi);
75 delete A_i;
76 p_mat_e_i->SetBlock(i, i, PtAP_i->EliminateCols(*ess_tdofs[i]));
77 PtAP_i->EliminateRows(*ess_tdofs[i]);
78 }
79 else
80 {
81 HypreParMatrix * Pj = (HypreParMatrix*)(&P->GetBlock(j,j));
82 A_r = new HypreParMatrix(trial_pfes[i]->GetComm(), trial_pfes[i]->GlobalVSize(),
83 trial_pfes[j]->GlobalVSize(), trial_pfes[i]->GetDofOffsets(),
84 trial_pfes[j]->GetDofOffsets(), &m_r->GetBlock(i,j));
85 PtAP_r = RAP(Pi,A_r,Pj);
86 delete A_r;
87 p_mat_e_r->SetBlock(i, j, PtAP_r->EliminateCols(*ess_tdofs[j]));
88 PtAP_r->EliminateRows(*ess_tdofs[i]);
89
90 A_i = new HypreParMatrix(trial_pfes[i]->GetComm(), trial_pfes[i]->GlobalVSize(),
91 trial_pfes[j]->GlobalVSize(), trial_pfes[i]->GetDofOffsets(),
92 trial_pfes[j]->GetDofOffsets(), &m_i->GetBlock(i,j));
93 PtAP_i = RAP(Pi,A_i,Pj);
94 delete A_i;
95 p_mat_e_i->SetBlock(i, j, PtAP_i->EliminateCols(*ess_tdofs[j]));
96 PtAP_i->EliminateRows(*ess_tdofs[i]);
97 }
98 p_mat_r->SetBlock(i, j, PtAP_r);
99 p_mat_i->SetBlock(i, j, PtAP_i);
100 }
101 }
102}
103
104
106{
109 P->owns_blocks = 0;
110 R->owns_blocks = 0;
111
112 for (int i = 0; i < nblocks; i++)
113 {
114 HypreParMatrix * P_ = trial_pfes[i]->Dof_TrueDof_Matrix();
115 P->SetBlock(i,i,P_);
116 const SparseMatrix * R_ = trial_pfes[i]->GetRestrictionMatrix();
117 R->SetBlock(i, i, const_cast<SparseMatrix*>(R_));
118 }
119}
120
122 &ess_tdof_list,
123 Vector &x,
125 Vector &X, Vector &B,
126 int copy_interior)
127{
128 FormSystemMatrix(ess_tdof_list, A);
129 if (static_cond)
130 {
131 static_cond->ReduceSystem(x, X, B, copy_interior);
132 }
133 else
134 {
135 int n = P->Width();
136 B.SetSize(2*n);
137 Vector B_r(B, 0, n);
138 Vector B_i(B, n, n);
139 P->MultTranspose(*y_r, B_r);
140 P->MultTranspose(*y_i, B_i);
141
142 int m = R->Height();
143 X.SetSize(2*m);
144
145 Vector X_r(X, 0, m);
146 Vector X_i(X, m, m);
147
148 Vector x_r(x, 0, x.Size()/2);
149 Vector x_i(x, x.Size()/2, x.Size()/2);
150
151 R->Mult(x_r, X_r);
152 R->Mult(x_i, X_i);
153
154 // eliminate tdof is RHS
155 // B_r -= Ae_r*X_r + Ae_i X_i
156 // B_i -= Ae_i*X_r + Ae_r X_i
157 Vector tmp(B_r.Size());
158 p_mat_e_r->Mult(X_r, tmp); B_r-=tmp;
159 p_mat_e_i->Mult(X_i, tmp); B_r+=tmp;
160
161 p_mat_e_i->Mult(X_r, tmp); B_i-=tmp;
162 p_mat_e_r->Mult(X_i, tmp); B_i-=tmp;
163
164 for (int j = 0; j < nblocks; j++)
165 {
166 if (!ess_tdofs[j]->Size()) { continue; }
167 for (int i = 0; i < ess_tdofs[j]->Size(); i++)
168 {
169 int tdof = (*ess_tdofs[j])[i];
170 int gdof = tdof + tdof_offsets[j];
171 B_r(gdof) = X_r(gdof); // diagonal policy is always one in parallel
172 B_i(gdof) = X_i(gdof); // diagonal policy is always one in parallel
173 }
174 }
175 if (!copy_interior)
176 {
177 X_r.SetSubVectorComplement(ess_tdof_list, 0.0);
178 X_i.SetSubVectorComplement(ess_tdof_list, 0.0);
179 }
180 }
181}
182
184 &ess_tdof_list,
186{
187 if (static_cond)
188 {
190 {
191 static_cond->SetEssentialTrueDofs(ess_tdof_list);
193 }
195 }
196 else
197 {
198 FillEssTdofLists(ess_tdof_list);
199 if (mat_r)
200 {
201 const int remove_zeros = 0;
202 Finalize(remove_zeros);
204 delete mat_r;
205 delete mat_i;
206 mat_r = nullptr;
207 mat_i = nullptr;
208 delete mat_e_r;
209 delete mat_e_i;
210 mat_e_r = nullptr;
211 mat_e_i = nullptr;
212 }
213 p_mat = new ComplexOperator(p_mat_r, p_mat_i, false, false);
214 A.Reset(p_mat, false);
215 }
216}
217
219 Vector &x)
220{
221 if (static_cond)
222 {
224 }
225 else
226 {
227 int n = P->Height();
228 int m = P->Width();
229 x.SetSize(2*n);
230
231 Vector x_r(x,0,n);
232 Vector x_i(x,n,n);
233 Vector X_r(const_cast<Vector&>(X), 0, m);
234 Vector X_i(const_cast<Vector&>(X), m, m);
235
236 P->Mult(X_r, x_r);
237 P->Mult(X_i, x_i);
238 }
239}
240
242{
244 delete p_mat_e_r;
245 delete p_mat_e_i;
246 p_mat_e_r = nullptr;
247 p_mat_e_i = nullptr;
248 delete p_mat_r;
249 delete p_mat_i;
250 p_mat_r = nullptr;
251 p_mat_i = nullptr;
252 delete p_mat;
253 p_mat = nullptr;
254 for (int i = 0; i < nblocks; i++)
255 {
256 delete ess_tdofs[i];
257 ess_tdofs[i] = new Array<int>();
258 }
259 delete P;
260 P = nullptr;
261 delete R;
262 R = nullptr;
263}
264
266{
267 delete p_mat_e_r;
268 delete p_mat_e_i;
269 p_mat_e_r = nullptr;
270 p_mat_e_i = nullptr;
271 delete p_mat_r;
272 delete p_mat_i;
273 p_mat_r = nullptr;
274 p_mat_i = nullptr;
275 delete p_mat;
276 p_mat = nullptr;
277 for (int i = 0; i < nblocks; i++)
278 {
279 delete ess_tdofs[i];
280 }
281 delete P;
282 delete R;
283}
284
285} // namespace mfem
286
287#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.
void ComputeSolution(const Vector &sc_sol, Vector &sol) const
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....
BlockMatrix * mat_e_r
Block Matrix used to store the eliminations from the b.c. Owned. .
BlockMatrix * mat_r
Block matrix to be associated with the real/imag Block bilinear form. Owned.
void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all integrators.
BlockVector * y_r
BlockVectors to be associated with the real/imag Block linear form.
ComplexBlockStaticCondensation * static_cond
Owned.
Mimic the action of a complex operator using two real operators.
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 * 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
void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Array< Array< int > * > ess_tdofs
ess_tdof list for each space
virtual ~ParComplexDPGWeakForm()
Destroys bilinear form.
Array< ParFiniteElementSpace * > trial_pfes
Trial FE spaces.
void FillEssTdofLists(const Array< int > &ess_tdof_list)
void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain integrators.
virtual void RecoverFEMSolution(const Vector &X, Vector &x)
void ParallelAssemble(BlockMatrix *mat_r, BlockMatrix *mat_i)
Returns the matrix assembled on the true dofs, i.e. P^t A P.
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)