1 // Copyright (c) 2010-2021, 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 "complex_fem.hpp"
13 #include "../general/forall.hpp"
14 
15 using namespace std;
16 
17 namespace mfem
18 {
19 
ComplexGridFunction(FiniteElementSpace * fes)20 ComplexGridFunction::ComplexGridFunction(FiniteElementSpace *fes)
21    : Vector(2*(fes->GetVSize()))
22 {
23    UseDevice(true);
24    this->Vector::operator=(0.0);
25 
26    gfr = new GridFunction();
27    gfr->MakeRef(fes, *this, 0);
28 
29    gfi = new GridFunction();
30    gfi->MakeRef(fes, *this, fes->GetVSize());
31 }
32 
33 void
Update()34 ComplexGridFunction::Update()
35 {
36    FiniteElementSpace *fes = gfr->FESpace();
37    const int vsize = fes->GetVSize();
38 
39    const Operator *T = fes->GetUpdateOperator();
40    if (T)
41    {
42       // Update the individual GridFunction objects. This will allocate new data
43       // arrays for each GridFunction.
44       gfr->Update();
45       gfi->Update();
46 
47       // Our data array now contains old data as well as being the wrong size so
48       // reallocate it.
49       UseDevice(true);
50       this->SetSize(2 * vsize);
51       this->Vector::operator=(0.0);
52 
53       // Create temporary vectors which point to the new data array
54       Vector gf_r; gf_r.MakeRef(*this, 0, vsize);
55       Vector gf_i; gf_i.MakeRef(*this, vsize, vsize);
56 
57       // Copy the updated GridFunctions into the new data array
58       gf_r = *gfr;
59       gf_i = *gfi;
60       gf_r.SyncAliasMemory(*this);
61       gf_i.SyncAliasMemory(*this);
62 
63       // Replace the individual data arrays with pointers into the new data
64       // array
65       gfr->MakeRef(*this, 0, vsize);
66       gfi->MakeRef(*this, vsize, vsize);
67    }
68    else
69    {
70       // The existing data will not be transferred to the new GridFunctions so
71       // delete it and allocate a new array
72       UseDevice(true);
73       this->SetSize(2 * vsize);
74       this->Vector::operator=(0.0);
75 
76       // Point the individual GridFunctions to the new data array
77       gfr->MakeRef(*this, 0, vsize);
78       gfi->MakeRef(*this, vsize, vsize);
79 
80       // These updates will only set the proper 'sequence' value within the
81       // individual GridFunction objects because their sizes are already correct
82       gfr->Update();
83       gfi->Update();
84    }
85 }
86 
87 void
ProjectCoefficient(Coefficient & real_coeff,Coefficient & imag_coeff)88 ComplexGridFunction::ProjectCoefficient(Coefficient &real_coeff,
89                                         Coefficient &imag_coeff)
90 {
91    gfr->SyncMemory(*this);
92    gfi->SyncMemory(*this);
93    gfr->ProjectCoefficient(real_coeff);
94    gfi->ProjectCoefficient(imag_coeff);
95    gfr->SyncAliasMemory(*this);
96    gfi->SyncAliasMemory(*this);
97 }
98 
99 void
ProjectCoefficient(VectorCoefficient & real_vcoeff,VectorCoefficient & imag_vcoeff)100 ComplexGridFunction::ProjectCoefficient(VectorCoefficient &real_vcoeff,
101                                         VectorCoefficient &imag_vcoeff)
102 {
103    gfr->SyncMemory(*this);
104    gfi->SyncMemory(*this);
105    gfr->ProjectCoefficient(real_vcoeff);
106    gfi->ProjectCoefficient(imag_vcoeff);
107    gfr->SyncAliasMemory(*this);
108    gfi->SyncAliasMemory(*this);
109 }
110 
111 void
ProjectBdrCoefficient(Coefficient & real_coeff,Coefficient & imag_coeff,Array<int> & attr)112 ComplexGridFunction::ProjectBdrCoefficient(Coefficient &real_coeff,
113                                            Coefficient &imag_coeff,
114                                            Array<int> &attr)
115 {
116    gfr->SyncMemory(*this);
117    gfi->SyncMemory(*this);
118    gfr->ProjectBdrCoefficient(real_coeff, attr);
119    gfi->ProjectBdrCoefficient(imag_coeff, attr);
120    gfr->SyncAliasMemory(*this);
121    gfi->SyncAliasMemory(*this);
122 }
123 
124 void
ProjectBdrCoefficientNormal(VectorCoefficient & real_vcoeff,VectorCoefficient & imag_vcoeff,Array<int> & attr)125 ComplexGridFunction::ProjectBdrCoefficientNormal(VectorCoefficient &real_vcoeff,
126                                                  VectorCoefficient &imag_vcoeff,
127                                                  Array<int> &attr)
128 {
129    gfr->SyncMemory(*this);
130    gfi->SyncMemory(*this);
131    gfr->ProjectBdrCoefficientNormal(real_vcoeff, attr);
132    gfi->ProjectBdrCoefficientNormal(imag_vcoeff, attr);
133    gfr->SyncAliasMemory(*this);
134    gfi->SyncAliasMemory(*this);
135 }
136 
137 void
ProjectBdrCoefficientTangent(VectorCoefficient & real_vcoeff,VectorCoefficient & imag_vcoeff,Array<int> & attr)138 ComplexGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient
139                                                   &real_vcoeff,
140                                                   VectorCoefficient
141                                                   &imag_vcoeff,
142                                                   Array<int> &attr)
143 {
144    gfr->SyncMemory(*this);
145    gfi->SyncMemory(*this);
146    gfr->ProjectBdrCoefficientTangent(real_vcoeff, attr);
147    gfi->ProjectBdrCoefficientTangent(imag_vcoeff, attr);
148    gfr->SyncAliasMemory(*this);
149    gfi->SyncAliasMemory(*this);
150 }
151 
152 
ComplexLinearForm(FiniteElementSpace * fes,ComplexOperator::Convention convention)153 ComplexLinearForm::ComplexLinearForm(FiniteElementSpace *fes,
154                                      ComplexOperator::Convention convention)
155    : Vector(2*(fes->GetVSize())),
156      conv(convention)
157 {
158    UseDevice(true);
159    this->Vector::operator=(0.0);
160 
161    lfr = new LinearForm();
162    lfr->MakeRef(fes, *this, 0);
163 
164    lfi = new LinearForm();
165    lfi->MakeRef(fes, *this, fes->GetVSize());
166 }
167 
ComplexLinearForm(FiniteElementSpace * fes,LinearForm * lf_r,LinearForm * lf_i,ComplexOperator::Convention convention)168 ComplexLinearForm::ComplexLinearForm(FiniteElementSpace *fes,
169                                      LinearForm *lf_r, LinearForm *lf_i,
170                                      ComplexOperator::Convention convention)
171    : Vector(2*(fes->GetVSize())),
172      conv(convention)
173 {
174    UseDevice(true);
175    this->Vector::operator=(0.0);
176 
177    lfr = new LinearForm(fes, lf_r);
178    lfi = new LinearForm(fes, lf_i);
179 
180    lfr->MakeRef(fes, *this, 0);
181    lfi->MakeRef(fes, *this, fes->GetVSize());
182 }
183 
~ComplexLinearForm()184 ComplexLinearForm::~ComplexLinearForm()
185 {
186    delete lfr;
187    delete lfi;
188 }
189 
190 void
AddDomainIntegrator(LinearFormIntegrator * lfi_real,LinearFormIntegrator * lfi_imag)191 ComplexLinearForm::AddDomainIntegrator(LinearFormIntegrator *lfi_real,
192                                        LinearFormIntegrator *lfi_imag)
193 {
194    if ( lfi_real ) { lfr->AddDomainIntegrator(lfi_real); }
195    if ( lfi_imag ) { lfi->AddDomainIntegrator(lfi_imag); }
196 }
197 
198 void
AddBoundaryIntegrator(LinearFormIntegrator * lfi_real,LinearFormIntegrator * lfi_imag)199 ComplexLinearForm::AddBoundaryIntegrator(LinearFormIntegrator *lfi_real,
200                                          LinearFormIntegrator *lfi_imag)
201 {
202    if ( lfi_real ) { lfr->AddBoundaryIntegrator(lfi_real); }
203    if ( lfi_imag ) { lfi->AddBoundaryIntegrator(lfi_imag); }
204 }
205 
206 void
AddBoundaryIntegrator(LinearFormIntegrator * lfi_real,LinearFormIntegrator * lfi_imag,Array<int> & bdr_attr_marker)207 ComplexLinearForm::AddBoundaryIntegrator(LinearFormIntegrator *lfi_real,
208                                          LinearFormIntegrator *lfi_imag,
209                                          Array<int> &bdr_attr_marker)
210 {
211    if ( lfi_real ) { lfr->AddBoundaryIntegrator(lfi_real, bdr_attr_marker); }
212    if ( lfi_imag ) { lfi->AddBoundaryIntegrator(lfi_imag, bdr_attr_marker); }
213 }
214 
215 void
AddBdrFaceIntegrator(LinearFormIntegrator * lfi_real,LinearFormIntegrator * lfi_imag)216 ComplexLinearForm::AddBdrFaceIntegrator(LinearFormIntegrator *lfi_real,
217                                         LinearFormIntegrator *lfi_imag)
218 {
219    if ( lfi_real ) { lfr->AddBdrFaceIntegrator(lfi_real); }
220    if ( lfi_imag ) { lfi->AddBdrFaceIntegrator(lfi_imag); }
221 }
222 
223 void
AddBdrFaceIntegrator(LinearFormIntegrator * lfi_real,LinearFormIntegrator * lfi_imag,Array<int> & bdr_attr_marker)224 ComplexLinearForm::AddBdrFaceIntegrator(LinearFormIntegrator *lfi_real,
225                                         LinearFormIntegrator *lfi_imag,
226                                         Array<int> &bdr_attr_marker)
227 {
228    if ( lfi_real ) { lfr->AddBdrFaceIntegrator(lfi_real, bdr_attr_marker); }
229    if ( lfi_imag ) { lfi->AddBdrFaceIntegrator(lfi_imag, bdr_attr_marker); }
230 }
231 
232 void
Update()233 ComplexLinearForm::Update()
234 {
235    FiniteElementSpace *fes = lfr->FESpace();
236    this->Update(fes);
237 }
238 
239 void
Update(FiniteElementSpace * fes)240 ComplexLinearForm::Update(FiniteElementSpace *fes)
241 {
242    UseDevice(true);
243    SetSize(2 * fes->GetVSize());
244    this->Vector::operator=(0.0);
245 
246    lfr->MakeRef(fes, *this, 0);
247    lfi->MakeRef(fes, *this, fes->GetVSize());
248 }
249 
250 void
Assemble()251 ComplexLinearForm::Assemble()
252 {
253    lfr->SyncMemory(*this);
254    lfi->SyncMemory(*this);
255    lfr->Assemble();
256    lfi->Assemble();
257    if (conv == ComplexOperator::BLOCK_SYMMETRIC) { *lfi *= -1.0; }
258    lfr->SyncAliasMemory(*this);
259    lfi->SyncAliasMemory(*this);
260 }
261 
262 complex<double>
operator ()(const ComplexGridFunction & gf) const263 ComplexLinearForm::operator()(const ComplexGridFunction &gf) const
264 {
265    double s = (conv == ComplexOperator::HERMITIAN) ? 1.0 : -1.0;
266    lfr->SyncMemory(*this);
267    lfi->SyncMemory(*this);
268    return complex<double>((*lfr)(gf.real()) - s * (*lfi)(gf.imag()),
269                           (*lfr)(gf.imag()) + s * (*lfi)(gf.real()));
270 }
271 
272 
RealInteg()273 bool SesquilinearForm::RealInteg()
274 {
275    int nint = blfr->GetFBFI()->Size() + blfr->GetDBFI()->Size() +
276               blfr->GetBBFI()->Size() + blfr->GetBFBFI()->Size();
277    return (nint != 0);
278 }
279 
ImagInteg()280 bool SesquilinearForm::ImagInteg()
281 {
282    int nint = blfi->GetFBFI()->Size() + blfi->GetDBFI()->Size() +
283               blfi->GetBBFI()->Size() + blfi->GetBFBFI()->Size();
284    return (nint != 0);
285 }
286 
SesquilinearForm(FiniteElementSpace * f,ComplexOperator::Convention convention)287 SesquilinearForm::SesquilinearForm(FiniteElementSpace *f,
288                                    ComplexOperator::Convention convention)
289    : conv(convention),
290      blfr(new BilinearForm(f)),
291      blfi(new BilinearForm(f))
292 {}
293 
SesquilinearForm(FiniteElementSpace * f,BilinearForm * bfr,BilinearForm * bfi,ComplexOperator::Convention convention)294 SesquilinearForm::SesquilinearForm(FiniteElementSpace *f,
295                                    BilinearForm *bfr, BilinearForm *bfi,
296                                    ComplexOperator::Convention convention)
297    : conv(convention),
298      blfr(new BilinearForm(f,bfr)),
299      blfi(new BilinearForm(f,bfi))
300 {}
301 
SetDiagonalPolicy(mfem::Matrix::DiagonalPolicy dpolicy)302 void SesquilinearForm::SetDiagonalPolicy(mfem::Matrix::DiagonalPolicy dpolicy)
303 {
304    diag_policy = dpolicy;
305 }
306 
~SesquilinearForm()307 SesquilinearForm::~SesquilinearForm()
308 {
309    delete blfr;
310    delete blfi;
311 }
312 
AddDomainIntegrator(BilinearFormIntegrator * bfi_real,BilinearFormIntegrator * bfi_imag)313 void SesquilinearForm::AddDomainIntegrator(BilinearFormIntegrator *bfi_real,
314                                            BilinearFormIntegrator *bfi_imag)
315 {
316    if (bfi_real) { blfr->AddDomainIntegrator(bfi_real); }
317    if (bfi_imag) { blfi->AddDomainIntegrator(bfi_imag); }
318 }
319 
320 void
AddBoundaryIntegrator(BilinearFormIntegrator * bfi_real,BilinearFormIntegrator * bfi_imag)321 SesquilinearForm::AddBoundaryIntegrator(BilinearFormIntegrator *bfi_real,
322                                         BilinearFormIntegrator *bfi_imag)
323 {
324    if (bfi_real) { blfr->AddBoundaryIntegrator(bfi_real); }
325    if (bfi_imag) { blfi->AddBoundaryIntegrator(bfi_imag); }
326 }
327 
328 void
AddBoundaryIntegrator(BilinearFormIntegrator * bfi_real,BilinearFormIntegrator * bfi_imag,Array<int> & bdr_marker)329 SesquilinearForm::AddBoundaryIntegrator(BilinearFormIntegrator *bfi_real,
330                                         BilinearFormIntegrator *bfi_imag,
331                                         Array<int> & bdr_marker)
332 {
333    if (bfi_real) { blfr->AddBoundaryIntegrator(bfi_real, bdr_marker); }
334    if (bfi_imag) { blfi->AddBoundaryIntegrator(bfi_imag, bdr_marker); }
335 }
336 
337 void
AddInteriorFaceIntegrator(BilinearFormIntegrator * bfi_real,BilinearFormIntegrator * bfi_imag)338 SesquilinearForm::AddInteriorFaceIntegrator(BilinearFormIntegrator *bfi_real,
339                                             BilinearFormIntegrator *bfi_imag)
340 {
341    if (bfi_real) { blfr->AddInteriorFaceIntegrator(bfi_real); }
342    if (bfi_imag) { blfi->AddInteriorFaceIntegrator(bfi_imag); }
343 }
344 
AddBdrFaceIntegrator(BilinearFormIntegrator * bfi_real,BilinearFormIntegrator * bfi_imag)345 void SesquilinearForm::AddBdrFaceIntegrator(BilinearFormIntegrator *bfi_real,
346                                             BilinearFormIntegrator *bfi_imag)
347 {
348    if (bfi_real) { blfr->AddBdrFaceIntegrator(bfi_real); }
349    if (bfi_imag) { blfi->AddBdrFaceIntegrator(bfi_imag); }
350 }
351 
AddBdrFaceIntegrator(BilinearFormIntegrator * bfi_real,BilinearFormIntegrator * bfi_imag,Array<int> & bdr_marker)352 void SesquilinearForm::AddBdrFaceIntegrator(BilinearFormIntegrator *bfi_real,
353                                             BilinearFormIntegrator *bfi_imag,
354                                             Array<int> &bdr_marker)
355 {
356    if (bfi_real) { blfr->AddBdrFaceIntegrator(bfi_real, bdr_marker); }
357    if (bfi_imag) { blfi->AddBdrFaceIntegrator(bfi_imag, bdr_marker); }
358 }
359 
360 void
Assemble(int skip_zeros)361 SesquilinearForm::Assemble(int skip_zeros)
362 {
363    blfr->Assemble(skip_zeros);
364    blfi->Assemble(skip_zeros);
365 }
366 
367 void
Finalize(int skip_zeros)368 SesquilinearForm::Finalize(int skip_zeros)
369 {
370    blfr->Finalize(skip_zeros);
371    blfi->Finalize(skip_zeros);
372 }
373 
374 ComplexSparseMatrix *
AssembleComplexSparseMatrix()375 SesquilinearForm::AssembleComplexSparseMatrix()
376 {
377    return new ComplexSparseMatrix(&blfr->SpMat(),
378                                   &blfi->SpMat(),
379                                   false, false, conv);
380 }
381 
382 void
FormLinearSystem(const Array<int> & ess_tdof_list,Vector & x,Vector & b,OperatorHandle & A,Vector & X,Vector & B,int ci)383 SesquilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list,
384                                    Vector &x, Vector &b,
385                                    OperatorHandle &A,
386                                    Vector &X, Vector &B,
387                                    int ci)
388 {
389    FiniteElementSpace *fes = blfr->FESpace();
390    const int vsize = fes->GetVSize();
391 
392    // Allocate temporary vector
393    Vector b_0;
394    b_0.UseDevice(true);
395    b_0.SetSize(vsize);
396    b_0 = 0.0;
397 
398    // Extract the real and imaginary parts of the input vectors
399    MFEM_ASSERT(x.Size() == 2 * vsize, "Input GridFunction of incorrect size!");
400    x.Read();
401    Vector x_r; x_r.MakeRef(x, 0, vsize);
402    Vector x_i; x_i.MakeRef(x, vsize, vsize);
403 
404    MFEM_ASSERT(b.Size() == 2 * vsize, "Input LinearForm of incorrect size!");
405    b.Read();
406    Vector b_r; b_r.MakeRef(b, 0, vsize);
407    Vector b_i; b_i.MakeRef(b, vsize, vsize);
408 
409    if (conv == ComplexOperator::BLOCK_SYMMETRIC) { b_i *= -1.0; }
410 
411    const int tvsize = fes->GetTrueVSize();
412    OperatorHandle A_r, A_i;
413 
414    X.UseDevice(true);
415    X.SetSize(2 * tvsize);
416    X = 0.0;
417 
418    B.UseDevice(true);
419    B.SetSize(2 * tvsize);
420    B = 0.0;
421 
422    Vector X_r; X_r.MakeRef(X, 0, tvsize);
423    Vector X_i; X_i.MakeRef(X, tvsize, tvsize);
424    Vector B_r; B_r.MakeRef(B, 0, tvsize);
425    Vector B_i; B_i.MakeRef(B, tvsize, tvsize);
426 
427    Vector X_0, B_0;
428 
429    if (RealInteg())
430    {
431       blfr->SetDiagonalPolicy(diag_policy);
432 
433       b_0 = b_r;
434       blfr->FormLinearSystem(ess_tdof_list, x_r, b_0, A_r, X_0, B_0, ci);
435       X_r = X_0; B_r = B_0;
436 
437       b_0 = b_i;
438       blfr->FormLinearSystem(ess_tdof_list, x_i, b_0, A_r, X_0, B_0, ci);
439       X_i = X_0; B_i = B_0;
440 
441       if (ImagInteg())
442       {
443          blfi->SetDiagonalPolicy(mfem::Matrix::DiagonalPolicy::DIAG_ZERO);
444 
445          b_0 = 0.0;
446          blfi->FormLinearSystem(ess_tdof_list, x_i, b_0, A_i, X_0, B_0, false);
447          B_r -= B_0;
448 
449          b_0 = 0.0;
450          blfi->FormLinearSystem(ess_tdof_list, x_r, b_0, A_i, X_0, B_0, false);
451          B_i += B_0;
452       }
453    }
454    else if (ImagInteg())
455    {
456       blfi->SetDiagonalPolicy(diag_policy);
457 
458       b_0 = b_i;
459       blfi->FormLinearSystem(ess_tdof_list, x_r, b_0, A_i, X_0, B_0, ci);
460       X_r = X_0; B_i = B_0;
461 
462       b_0 = b_r; b_0 *= -1.0;
463       blfi->FormLinearSystem(ess_tdof_list, x_i, b_0, A_i, X_0, B_0, ci);
464       X_i = X_0; B_r = B_0; B_r *= -1.0;
465    }
466    else
467    {
468       MFEM_ABORT("Real and Imaginary part of the Sesquilinear form are empty");
469    }
470 
471    if (RealInteg() && ImagInteg())
472    {
473       // Modify RHS and offdiagonal blocks (imaginary parts of the matrix) to
474       // conform with standard essential BC treatment
475       if (A_i.Is<ConstrainedOperator>())
476       {
477          const int n = ess_tdof_list.Size();
478          auto d_B_r = B_r.Write();
479          auto d_B_i = B_i.Write();
480          auto d_X_r = X_r.Read();
481          auto d_X_i = X_i.Read();
482          auto d_idx = ess_tdof_list.Read();
483          MFEM_FORALL(i, n,
484          {
485             const int j = d_idx[i];
486             d_B_r[j] = d_X_r[j];
487             d_B_i[j] = d_X_i[j];
488          });
489          A_i.As<ConstrainedOperator>()->SetDiagonalPolicy
490          (mfem::Operator::DiagonalPolicy::DIAG_ZERO);
491       }
492    }
493 
494    if (conv == ComplexOperator::BLOCK_SYMMETRIC)
495    {
496       B_i *= -1.0;
497       b_i *= -1.0;
498    }
499 
500    x_r.SyncAliasMemory(x);
501    x_i.SyncAliasMemory(x);
502    b_r.SyncAliasMemory(b);
503    b_i.SyncAliasMemory(b);
504 
505    X_r.SyncAliasMemory(X);
506    X_i.SyncAliasMemory(X);
507    B_r.SyncAliasMemory(B);
508    B_i.SyncAliasMemory(B);
509 
510    // A = A_r + i A_i
511    A.Clear();
512    if ( A_r.Type() == Operator::MFEM_SPARSEMAT ||
513         A_i.Type() == Operator::MFEM_SPARSEMAT )
514    {
515       ComplexSparseMatrix * A_sp =
516          new ComplexSparseMatrix(A_r.As<SparseMatrix>(),
517                                  A_i.As<SparseMatrix>(),
518                                  A_r.OwnsOperator(),
519                                  A_i.OwnsOperator(),
520                                  conv);
521       A.Reset<ComplexSparseMatrix>(A_sp, true);
522    }
523    else
524    {
525       ComplexOperator * A_op =
526          new ComplexOperator(A_r.Ptr(),
527                              A_i.Ptr(),
528                              A_r.OwnsOperator(),
529                              A_i.OwnsOperator(),
530                              conv);
531       A.Reset<ComplexOperator>(A_op, true);
532    }
533    A_r.SetOperatorOwner(false);
534    A_i.SetOperatorOwner(false);
535 }
536 
537 void
FormSystemMatrix(const Array<int> & ess_tdof_list,OperatorHandle & A)538 SesquilinearForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
539                                    OperatorHandle &A)
540 
541 {
542    OperatorHandle A_r, A_i;
543    if (RealInteg())
544    {
545       blfr->SetDiagonalPolicy(diag_policy);
546       blfr->FormSystemMatrix(ess_tdof_list, A_r);
547    }
548    if (ImagInteg())
549    {
550       blfi->SetDiagonalPolicy(RealInteg() ?
551                               mfem::Matrix::DiagonalPolicy::DIAG_ZERO :
552                               diag_policy);
553       blfi->FormSystemMatrix(ess_tdof_list, A_i);
554    }
555    if (!RealInteg() && !ImagInteg())
556    {
557       MFEM_ABORT("Both Real and Imaginary part of the Sesquilinear form are empty");
558    }
559 
560    if (RealInteg() && ImagInteg())
561    {
562       // Modify offdiagonal blocks (imaginary parts of the matrix) to conform
563       // with standard essential BC treatment
564       if (A_i.Is<ConstrainedOperator>())
565       {
566          A_i.As<ConstrainedOperator>()->SetDiagonalPolicy
567          (mfem::Operator::DiagonalPolicy::DIAG_ZERO);
568       }
569    }
570 
571    // A = A_r + i A_i
572    A.Clear();
573    if ( A_r.Type() == Operator::MFEM_SPARSEMAT ||
574         A_i.Type() == Operator::MFEM_SPARSEMAT )
575    {
576       ComplexSparseMatrix * A_sp =
577          new ComplexSparseMatrix(A_r.As<SparseMatrix>(),
578                                  A_i.As<SparseMatrix>(),
579                                  A_r.OwnsOperator(),
580                                  A_i.OwnsOperator(),
581                                  conv);
582       A.Reset<ComplexSparseMatrix>(A_sp, true);
583    }
584    else
585    {
586       ComplexOperator * A_op =
587          new ComplexOperator(A_r.Ptr(),
588                              A_i.Ptr(),
589                              A_r.OwnsOperator(),
590                              A_i.OwnsOperator(),
591                              conv);
592       A.Reset<ComplexOperator>(A_op, true);
593    }
594    A_r.SetOperatorOwner(false);
595    A_i.SetOperatorOwner(false);
596 }
597 
598 void
RecoverFEMSolution(const Vector & X,const Vector & b,Vector & x)599 SesquilinearForm::RecoverFEMSolution(const Vector &X, const Vector &b,
600                                      Vector &x)
601 {
602    FiniteElementSpace *fes = blfr->FESpace();
603 
604    const SparseMatrix *P = fes->GetConformingProlongation();
605    if (!P)
606    {
607       x = X;
608       return;
609    }
610 
611    const int vsize  = fes->GetVSize();
612    const int tvsize = X.Size() / 2;
613 
614    X.Read();
615    Vector X_r; X_r.MakeRef(const_cast<Vector&>(X), 0, tvsize);
616    Vector X_i; X_i.MakeRef(const_cast<Vector&>(X), tvsize, tvsize);
617 
618    x.Write();
619    Vector x_r; x_r.MakeRef(x, 0, vsize);
620    Vector x_i; x_i.MakeRef(x, vsize, vsize);
621 
622    // Apply conforming prolongation
623    P->Mult(X_r, x_r);
624    P->Mult(X_i, x_i);
625 
626    x_r.SyncAliasMemory(x);
627    x_i.SyncAliasMemory(x);
628 }
629 
630 void
Update(FiniteElementSpace * nfes)631 SesquilinearForm::Update(FiniteElementSpace *nfes)
632 {
633    if ( blfr ) { blfr->Update(nfes); }
634    if ( blfi ) { blfi->Update(nfes); }
635 }
636 
637 
638 #ifdef MFEM_USE_MPI
639 
ParComplexGridFunction(ParFiniteElementSpace * pfes)640 ParComplexGridFunction::ParComplexGridFunction(ParFiniteElementSpace *pfes)
641    : Vector(2*(pfes->GetVSize()))
642 {
643    UseDevice(true);
644    this->Vector::operator=(0.0);
645 
646    pgfr = new ParGridFunction();
647    pgfr->MakeRef(pfes, *this, 0);
648 
649    pgfi = new ParGridFunction();
650    pgfi->MakeRef(pfes, *this, pfes->GetVSize());
651 }
652 
653 void
Update()654 ParComplexGridFunction::Update()
655 {
656    ParFiniteElementSpace *pfes = pgfr->ParFESpace();
657    const int vsize = pfes->GetVSize();
658 
659    const Operator *T = pfes->GetUpdateOperator();
660    if (T)
661    {
662       // Update the individual GridFunction objects. This will allocate new data
663       // arrays for each GridFunction.
664       pgfr->Update();
665       pgfi->Update();
666 
667       // Our data array now contains old data as well as being the wrong size so
668       // reallocate it.
669       UseDevice(true);
670       this->SetSize(2 * vsize);
671       this->Vector::operator=(0.0);
672 
673       // Create temporary vectors which point to the new data array
674       Vector gf_r; gf_r.MakeRef(*this, 0, vsize);
675       Vector gf_i; gf_i.MakeRef(*this, vsize, vsize);
676 
677       // Copy the updated GridFunctions into the new data array
678       gf_r = *pgfr; gf_r.SyncAliasMemory(*this);
679       gf_i = *pgfi; gf_i.SyncAliasMemory(*this);
680 
681       // Replace the individual data arrays with pointers into the new data
682       // array
683       pgfr->MakeRef(*this, 0, vsize);
684       pgfi->MakeRef(*this, vsize, vsize);
685    }
686    else
687    {
688       // The existing data will not be transferred to the new GridFunctions so
689       // delete it and allocate a new array
690       UseDevice(true);
691       this->SetSize(2 * vsize);
692       this->Vector::operator=(0.0);
693 
694       // Point the individual GridFunctions to the new data array
695       pgfr->MakeRef(*this, 0, vsize);
696       pgfi->MakeRef(*this, vsize, vsize);
697 
698       // These updates will only set the proper 'sequence' value within the
699       // individual GridFunction objects because their sizes are already correct
700       pgfr->Update();
701       pgfi->Update();
702    }
703 }
704 
705 void
ProjectCoefficient(Coefficient & real_coeff,Coefficient & imag_coeff)706 ParComplexGridFunction::ProjectCoefficient(Coefficient &real_coeff,
707                                            Coefficient &imag_coeff)
708 {
709    pgfr->SyncMemory(*this);
710    pgfi->SyncMemory(*this);
711    pgfr->ProjectCoefficient(real_coeff);
712    pgfi->ProjectCoefficient(imag_coeff);
713    pgfr->SyncAliasMemory(*this);
714    pgfi->SyncAliasMemory(*this);
715 }
716 
717 void
ProjectCoefficient(VectorCoefficient & real_vcoeff,VectorCoefficient & imag_vcoeff)718 ParComplexGridFunction::ProjectCoefficient(VectorCoefficient &real_vcoeff,
719                                            VectorCoefficient &imag_vcoeff)
720 {
721    pgfr->SyncMemory(*this);
722    pgfi->SyncMemory(*this);
723    pgfr->ProjectCoefficient(real_vcoeff);
724    pgfi->ProjectCoefficient(imag_vcoeff);
725    pgfr->SyncAliasMemory(*this);
726    pgfi->SyncAliasMemory(*this);
727 }
728 
729 void
ProjectBdrCoefficient(Coefficient & real_coeff,Coefficient & imag_coeff,Array<int> & attr)730 ParComplexGridFunction::ProjectBdrCoefficient(Coefficient &real_coeff,
731                                               Coefficient &imag_coeff,
732                                               Array<int> &attr)
733 {
734    pgfr->SyncMemory(*this);
735    pgfi->SyncMemory(*this);
736    pgfr->ProjectBdrCoefficient(real_coeff, attr);
737    pgfi->ProjectBdrCoefficient(imag_coeff, attr);
738    pgfr->SyncAliasMemory(*this);
739    pgfi->SyncAliasMemory(*this);
740 }
741 
742 void
ProjectBdrCoefficientNormal(VectorCoefficient & real_vcoeff,VectorCoefficient & imag_vcoeff,Array<int> & attr)743 ParComplexGridFunction::ProjectBdrCoefficientNormal(VectorCoefficient
744                                                     &real_vcoeff,
745                                                     VectorCoefficient
746                                                     &imag_vcoeff,
747                                                     Array<int> &attr)
748 {
749    pgfr->SyncMemory(*this);
750    pgfi->SyncMemory(*this);
751    pgfr->ProjectBdrCoefficientNormal(real_vcoeff, attr);
752    pgfi->ProjectBdrCoefficientNormal(imag_vcoeff, attr);
753    pgfr->SyncAliasMemory(*this);
754    pgfi->SyncAliasMemory(*this);
755 }
756 
757 void
ProjectBdrCoefficientTangent(VectorCoefficient & real_vcoeff,VectorCoefficient & imag_vcoeff,Array<int> & attr)758 ParComplexGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient
759                                                      &real_vcoeff,
760                                                      VectorCoefficient
761                                                      &imag_vcoeff,
762                                                      Array<int> &attr)
763 {
764    pgfr->SyncMemory(*this);
765    pgfi->SyncMemory(*this);
766    pgfr->ProjectBdrCoefficientTangent(real_vcoeff, attr);
767    pgfi->ProjectBdrCoefficientTangent(imag_vcoeff, attr);
768    pgfr->SyncAliasMemory(*this);
769    pgfi->SyncAliasMemory(*this);
770 }
771 
772 void
Distribute(const Vector * tv)773 ParComplexGridFunction::Distribute(const Vector *tv)
774 {
775    ParFiniteElementSpace *pfes = pgfr->ParFESpace();
776    const int tvsize = pfes->GetTrueVSize();
777 
778    tv->Read();
779    Vector tvr; tvr.MakeRef(const_cast<Vector&>(*tv), 0, tvsize);
780    Vector tvi; tvi.MakeRef(const_cast<Vector&>(*tv), tvsize, tvsize);
781 
782    pgfr->SyncMemory(*this);
783    pgfi->SyncMemory(*this);
784    pgfr->Distribute(tvr);
785    pgfi->Distribute(tvi);
786    pgfr->SyncAliasMemory(*this);
787    pgfi->SyncAliasMemory(*this);
788 }
789 
790 void
ParallelProject(Vector & tv) const791 ParComplexGridFunction::ParallelProject(Vector &tv) const
792 {
793    ParFiniteElementSpace *pfes = pgfr->ParFESpace();
794    const int tvsize = pfes->GetTrueVSize();
795 
796    tv.Write();
797    Vector tvr; tvr.MakeRef(tv, 0, tvsize);
798    Vector tvi; tvi.MakeRef(tv, tvsize, tvsize);
799 
800    pgfr->SyncMemory(*this);
801    pgfi->SyncMemory(*this);
802    pgfr->ParallelProject(tvr);
803    pgfi->ParallelProject(tvi);
804    pgfr->SyncAliasMemory(*this);
805    pgfi->SyncAliasMemory(*this);
806 
807    tvr.SyncAliasMemory(tv);
808    tvi.SyncAliasMemory(tv);
809 }
810 
811 
ParComplexLinearForm(ParFiniteElementSpace * pfes,ComplexOperator::Convention convention)812 ParComplexLinearForm::ParComplexLinearForm(ParFiniteElementSpace *pfes,
813                                            ComplexOperator::Convention
814                                            convention)
815    : Vector(2*(pfes->GetVSize())),
816      conv(convention)
817 {
818    UseDevice(true);
819    this->Vector::operator=(0.0);
820 
821    plfr = new ParLinearForm();
822    plfr->MakeRef(pfes, *this, 0);
823 
824    plfi = new ParLinearForm();
825    plfi->MakeRef(pfes, *this, pfes->GetVSize());
826 
827    HYPRE_BigInt *tdof_offsets_fes = pfes->GetTrueDofOffsets();
828 
829    int n = (HYPRE_AssumedPartitionCheck()) ? 2 : pfes->GetNRanks();
830    tdof_offsets = new HYPRE_BigInt[n+1];
831 
832    for (int i = 0; i <= n; i++)
833    {
834       tdof_offsets[i] = 2 * tdof_offsets_fes[i];
835    }
836 }
837 
838 
ParComplexLinearForm(ParFiniteElementSpace * pfes,ParLinearForm * plf_r,ParLinearForm * plf_i,ComplexOperator::Convention convention)839 ParComplexLinearForm::ParComplexLinearForm(ParFiniteElementSpace *pfes,
840                                            ParLinearForm *plf_r,
841                                            ParLinearForm *plf_i,
842                                            ComplexOperator::Convention
843                                            convention)
844    : Vector(2*(pfes->GetVSize())),
845      conv(convention)
846 {
847    UseDevice(true);
848    this->Vector::operator=(0.0);
849 
850    plfr = new ParLinearForm(pfes, plf_r);
851    plfi = new ParLinearForm(pfes, plf_i);
852 
853    plfr->MakeRef(pfes, *this, 0);
854    plfi->MakeRef(pfes, *this, pfes->GetVSize());
855 
856    HYPRE_BigInt *tdof_offsets_fes = pfes->GetTrueDofOffsets();
857 
858    int n = (HYPRE_AssumedPartitionCheck()) ? 2 : pfes->GetNRanks();
859    tdof_offsets = new HYPRE_BigInt[n+1];
860 
861    for (int i = 0; i <= n; i++)
862    {
863       tdof_offsets[i] = 2 * tdof_offsets_fes[i];
864    }
865 }
866 
~ParComplexLinearForm()867 ParComplexLinearForm::~ParComplexLinearForm()
868 {
869    delete plfr;
870    delete plfi;
871    delete [] tdof_offsets;
872 }
873 
874 void
AddDomainIntegrator(LinearFormIntegrator * lfi_real,LinearFormIntegrator * lfi_imag)875 ParComplexLinearForm::AddDomainIntegrator(LinearFormIntegrator *lfi_real,
876                                           LinearFormIntegrator *lfi_imag)
877 {
878    if ( lfi_real ) { plfr->AddDomainIntegrator(lfi_real); }
879    if ( lfi_imag ) { plfi->AddDomainIntegrator(lfi_imag); }
880 }
881 
882 void
AddBoundaryIntegrator(LinearFormIntegrator * lfi_real,LinearFormIntegrator * lfi_imag)883 ParComplexLinearForm::AddBoundaryIntegrator(LinearFormIntegrator *lfi_real,
884                                             LinearFormIntegrator *lfi_imag)
885 {
886    if ( lfi_real ) { plfr->AddBoundaryIntegrator(lfi_real); }
887    if ( lfi_imag ) { plfi->AddBoundaryIntegrator(lfi_imag); }
888 }
889 
890 void
AddBoundaryIntegrator(LinearFormIntegrator * lfi_real,LinearFormIntegrator * lfi_imag,Array<int> & bdr_attr_marker)891 ParComplexLinearForm::AddBoundaryIntegrator(LinearFormIntegrator *lfi_real,
892                                             LinearFormIntegrator *lfi_imag,
893                                             Array<int> &bdr_attr_marker)
894 {
895    if ( lfi_real ) { plfr->AddBoundaryIntegrator(lfi_real, bdr_attr_marker); }
896    if ( lfi_imag ) { plfi->AddBoundaryIntegrator(lfi_imag, bdr_attr_marker); }
897 }
898 
899 void
AddBdrFaceIntegrator(LinearFormIntegrator * lfi_real,LinearFormIntegrator * lfi_imag)900 ParComplexLinearForm::AddBdrFaceIntegrator(LinearFormIntegrator *lfi_real,
901                                            LinearFormIntegrator *lfi_imag)
902 {
903    if ( lfi_real ) { plfr->AddBdrFaceIntegrator(lfi_real); }
904    if ( lfi_imag ) { plfi->AddBdrFaceIntegrator(lfi_imag); }
905 }
906 
907 void
AddBdrFaceIntegrator(LinearFormIntegrator * lfi_real,LinearFormIntegrator * lfi_imag,Array<int> & bdr_attr_marker)908 ParComplexLinearForm::AddBdrFaceIntegrator(LinearFormIntegrator *lfi_real,
909                                            LinearFormIntegrator *lfi_imag,
910                                            Array<int> &bdr_attr_marker)
911 {
912    if ( lfi_real ) { plfr->AddBdrFaceIntegrator(lfi_real, bdr_attr_marker); }
913    if ( lfi_imag ) { plfi->AddBdrFaceIntegrator(lfi_imag, bdr_attr_marker); }
914 }
915 
916 void
Update(ParFiniteElementSpace * pf)917 ParComplexLinearForm::Update(ParFiniteElementSpace *pf)
918 {
919    ParFiniteElementSpace *pfes = (pf != NULL) ? pf : plfr->ParFESpace();
920 
921    UseDevice(true);
922    SetSize(2 * pfes->GetVSize());
923    this->Vector::operator=(0.0);
924 
925    plfr->MakeRef(pfes, *this, 0);
926    plfi->MakeRef(pfes, *this, pfes->GetVSize());
927 }
928 
929 void
Assemble()930 ParComplexLinearForm::Assemble()
931 {
932    plfr->SyncMemory(*this);
933    plfi->SyncMemory(*this);
934    plfr->Assemble();
935    plfi->Assemble();
936    if (conv == ComplexOperator::BLOCK_SYMMETRIC) { *plfi *= -1.0; }
937    plfr->SyncAliasMemory(*this);
938    plfi->SyncAliasMemory(*this);
939 }
940 
941 void
ParallelAssemble(Vector & tv)942 ParComplexLinearForm::ParallelAssemble(Vector &tv)
943 {
944    const int tvsize = plfr->ParFESpace()->GetTrueVSize();
945 
946    tv.Write();
947    Vector tvr; tvr.MakeRef(tv, 0, tvsize);
948    Vector tvi; tvi.MakeRef(tv, tvsize, tvsize);
949 
950    plfr->SyncMemory(*this);
951    plfi->SyncMemory(*this);
952    plfr->ParallelAssemble(tvr);
953    plfi->ParallelAssemble(tvi);
954    plfr->SyncAliasMemory(*this);
955    plfi->SyncAliasMemory(*this);
956 
957    tvr.SyncAliasMemory(tv);
958    tvi.SyncAliasMemory(tv);
959 }
960 
961 HypreParVector *
ParallelAssemble()962 ParComplexLinearForm::ParallelAssemble()
963 {
964    const ParFiniteElementSpace *pfes = plfr->ParFESpace();
965    const int tvsize = pfes->GetTrueVSize();
966 
967    HypreParVector *tv = new HypreParVector(pfes->GetComm(),
968                                            2*(pfes->GlobalTrueVSize()),
969                                            tdof_offsets);
970 
971    tv->Write();
972    Vector tvr; tvr.MakeRef(*tv, 0, tvsize);
973    Vector tvi; tvi.MakeRef(*tv, tvsize, tvsize);
974 
975    plfr->SyncMemory(*this);
976    plfi->SyncMemory(*this);
977    plfr->ParallelAssemble(tvr);
978    plfi->ParallelAssemble(tvi);
979    plfr->SyncAliasMemory(*this);
980    plfi->SyncAliasMemory(*this);
981 
982    tvr.SyncAliasMemory(*tv);
983    tvi.SyncAliasMemory(*tv);
984 
985    return tv;
986 }
987 
988 complex<double>
operator ()(const ParComplexGridFunction & gf) const989 ParComplexLinearForm::operator()(const ParComplexGridFunction &gf) const
990 {
991    plfr->SyncMemory(*this);
992    plfi->SyncMemory(*this);
993    double s = (conv == ComplexOperator::HERMITIAN) ? 1.0 : -1.0;
994    return complex<double>((*plfr)(gf.real()) - s * (*plfi)(gf.imag()),
995                           (*plfr)(gf.imag()) + s * (*plfi)(gf.real()));
996 }
997 
998 
RealInteg()999 bool ParSesquilinearForm::RealInteg()
1000 {
1001    int nint = pblfr->GetFBFI()->Size() + pblfr->GetDBFI()->Size() +
1002               pblfr->GetBBFI()->Size() + pblfr->GetBFBFI()->Size();
1003    return (nint != 0);
1004 }
1005 
ImagInteg()1006 bool ParSesquilinearForm::ImagInteg()
1007 {
1008    int nint = pblfi->GetFBFI()->Size() + pblfi->GetDBFI()->Size() +
1009               pblfi->GetBBFI()->Size() + pblfi->GetBFBFI()->Size();
1010    return (nint != 0);
1011 }
1012 
ParSesquilinearForm(ParFiniteElementSpace * pf,ComplexOperator::Convention convention)1013 ParSesquilinearForm::ParSesquilinearForm(ParFiniteElementSpace *pf,
1014                                          ComplexOperator::Convention
1015                                          convention)
1016    : conv(convention),
1017      pblfr(new ParBilinearForm(pf)),
1018      pblfi(new ParBilinearForm(pf))
1019 {}
1020 
ParSesquilinearForm(ParFiniteElementSpace * pf,ParBilinearForm * pbfr,ParBilinearForm * pbfi,ComplexOperator::Convention convention)1021 ParSesquilinearForm::ParSesquilinearForm(ParFiniteElementSpace *pf,
1022                                          ParBilinearForm *pbfr,
1023                                          ParBilinearForm *pbfi,
1024                                          ComplexOperator::Convention convention)
1025    : conv(convention),
1026      pblfr(new ParBilinearForm(pf,pbfr)),
1027      pblfi(new ParBilinearForm(pf,pbfi))
1028 {}
1029 
~ParSesquilinearForm()1030 ParSesquilinearForm::~ParSesquilinearForm()
1031 {
1032    delete pblfr;
1033    delete pblfi;
1034 }
1035 
AddDomainIntegrator(BilinearFormIntegrator * bfi_real,BilinearFormIntegrator * bfi_imag)1036 void ParSesquilinearForm::AddDomainIntegrator(BilinearFormIntegrator *bfi_real,
1037                                               BilinearFormIntegrator *bfi_imag)
1038 {
1039    if (bfi_real) { pblfr->AddDomainIntegrator(bfi_real); }
1040    if (bfi_imag) { pblfi->AddDomainIntegrator(bfi_imag); }
1041 }
1042 
1043 void
AddBoundaryIntegrator(BilinearFormIntegrator * bfi_real,BilinearFormIntegrator * bfi_imag)1044 ParSesquilinearForm::AddBoundaryIntegrator(BilinearFormIntegrator *bfi_real,
1045                                            BilinearFormIntegrator *bfi_imag)
1046 {
1047    if (bfi_real) { pblfr->AddBoundaryIntegrator(bfi_real); }
1048    if (bfi_imag) { pblfi->AddBoundaryIntegrator(bfi_imag); }
1049 }
1050 
1051 void
AddBoundaryIntegrator(BilinearFormIntegrator * bfi_real,BilinearFormIntegrator * bfi_imag,Array<int> & bdr_marker)1052 ParSesquilinearForm::AddBoundaryIntegrator(BilinearFormIntegrator *bfi_real,
1053                                            BilinearFormIntegrator *bfi_imag,
1054                                            Array<int> & bdr_marker)
1055 {
1056    if (bfi_real) { pblfr->AddBoundaryIntegrator(bfi_real, bdr_marker); }
1057    if (bfi_imag) { pblfi->AddBoundaryIntegrator(bfi_imag, bdr_marker); }
1058 }
1059 
1060 void
AddInteriorFaceIntegrator(BilinearFormIntegrator * bfi_real,BilinearFormIntegrator * bfi_imag)1061 ParSesquilinearForm::AddInteriorFaceIntegrator(BilinearFormIntegrator *bfi_real,
1062                                                BilinearFormIntegrator *bfi_imag)
1063 {
1064    if (bfi_real) { pblfr->AddInteriorFaceIntegrator(bfi_real); }
1065    if (bfi_imag) { pblfi->AddInteriorFaceIntegrator(bfi_imag); }
1066 }
1067 
1068 void
AddBdrFaceIntegrator(BilinearFormIntegrator * bfi_real,BilinearFormIntegrator * bfi_imag)1069 ParSesquilinearForm::AddBdrFaceIntegrator(BilinearFormIntegrator *bfi_real,
1070                                           BilinearFormIntegrator *bfi_imag)
1071 {
1072    if (bfi_real) { pblfr->AddBdrFaceIntegrator(bfi_real); }
1073    if (bfi_imag) { pblfi->AddBdrFaceIntegrator(bfi_imag); }
1074 }
1075 
1076 void
AddBdrFaceIntegrator(BilinearFormIntegrator * bfi_real,BilinearFormIntegrator * bfi_imag,Array<int> & bdr_marker)1077 ParSesquilinearForm::AddBdrFaceIntegrator(BilinearFormIntegrator *bfi_real,
1078                                           BilinearFormIntegrator *bfi_imag,
1079                                           Array<int> &bdr_marker)
1080 {
1081    if (bfi_real) { pblfr->AddBdrFaceIntegrator(bfi_real, bdr_marker); }
1082    if (bfi_imag) { pblfi->AddBdrFaceIntegrator(bfi_imag, bdr_marker); }
1083 }
1084 
1085 void
Assemble(int skip_zeros)1086 ParSesquilinearForm::Assemble(int skip_zeros)
1087 {
1088    pblfr->Assemble(skip_zeros);
1089    pblfi->Assemble(skip_zeros);
1090 }
1091 
1092 void
Finalize(int skip_zeros)1093 ParSesquilinearForm::Finalize(int skip_zeros)
1094 {
1095    pblfr->Finalize(skip_zeros);
1096    pblfi->Finalize(skip_zeros);
1097 }
1098 
1099 ComplexHypreParMatrix *
ParallelAssemble()1100 ParSesquilinearForm::ParallelAssemble()
1101 {
1102    return new ComplexHypreParMatrix(pblfr->ParallelAssemble(),
1103                                     pblfi->ParallelAssemble(),
1104                                     true, true, conv);
1105 }
1106 
1107 void
FormLinearSystem(const Array<int> & ess_tdof_list,Vector & x,Vector & b,OperatorHandle & A,Vector & X,Vector & B,int ci)1108 ParSesquilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list,
1109                                       Vector &x, Vector &b,
1110                                       OperatorHandle &A,
1111                                       Vector &X, Vector &B,
1112                                       int ci)
1113 {
1114    ParFiniteElementSpace *pfes = pblfr->ParFESpace();
1115    const int vsize = pfes->GetVSize();
1116 
1117    // Allocate temporary vector
1118    Vector b_0;
1119    b_0.UseDevice(true);
1120    b_0.SetSize(vsize);
1121    b_0 = 0.0;
1122 
1123    // Extract the real and imaginary parts of the input vectors
1124    MFEM_ASSERT(x.Size() == 2 * vsize, "Input GridFunction of incorrect size!");
1125    x.Read();
1126    Vector x_r; x_r.MakeRef(x, 0, vsize);
1127    Vector x_i; x_i.MakeRef(x, vsize, vsize);
1128 
1129    MFEM_ASSERT(b.Size() == 2 * vsize, "Input LinearForm of incorrect size!");
1130    b.Read();
1131    Vector b_r; b_r.MakeRef(b, 0, vsize);
1132    Vector b_i; b_i.MakeRef(b, vsize, vsize);
1133 
1134    if (conv == ComplexOperator::BLOCK_SYMMETRIC) { b_i *= -1.0; }
1135 
1136    const int tvsize = pfes->GetTrueVSize();
1137    OperatorHandle A_r, A_i;
1138 
1139    X.UseDevice(true);
1140    X.SetSize(2 * tvsize);
1141    X = 0.0;
1142 
1143    B.UseDevice(true);
1144    B.SetSize(2 * tvsize);
1145    B = 0.0;
1146 
1147    Vector X_r; X_r.MakeRef(X, 0, tvsize);
1148    Vector X_i; X_i.MakeRef(X, tvsize, tvsize);
1149    Vector B_r; B_r.MakeRef(B, 0, tvsize);
1150    Vector B_i; B_i.MakeRef(B, tvsize, tvsize);
1151 
1152    Vector X_0, B_0;
1153 
1154    if (RealInteg())
1155    {
1156       b_0 = b_r;
1157       pblfr->FormLinearSystem(ess_tdof_list, x_r, b_0, A_r, X_0, B_0, ci);
1158       X_r = X_0; B_r = B_0;
1159 
1160       b_0 = b_i;
1161       pblfr->FormLinearSystem(ess_tdof_list, x_i, b_0, A_r, X_0, B_0, ci);
1162       X_i = X_0; B_i = B_0;
1163 
1164       if (ImagInteg())
1165       {
1166          b_0 = 0.0;
1167          pblfi->FormLinearSystem(ess_tdof_list, x_i, b_0, A_i, X_0, B_0, false);
1168          B_r -= B_0;
1169 
1170          b_0 = 0.0;
1171          pblfi->FormLinearSystem(ess_tdof_list, x_r, b_0, A_i, X_0, B_0, false);
1172          B_i += B_0;
1173       }
1174    }
1175    else if (ImagInteg())
1176    {
1177       b_0 = b_i;
1178       pblfi->FormLinearSystem(ess_tdof_list, x_r, b_0, A_i, X_0, B_0, ci);
1179       X_r = X_0; B_i = B_0;
1180 
1181       b_0 = b_r; b_0 *= -1.0;
1182       pblfi->FormLinearSystem(ess_tdof_list, x_i, b_0, A_i, X_0, B_0, ci);
1183       X_i = X_0; B_r = B_0; B_r *= -1.0;
1184    }
1185    else
1186    {
1187       MFEM_ABORT("Real and Imaginary part of the Sesquilinear form are empty");
1188    }
1189 
1190    if (RealInteg() && ImagInteg())
1191    {
1192       // Modify RHS to conform with standard essential BC treatment
1193       const int n = ess_tdof_list.Size();
1194       auto d_B_r = B_r.Write();
1195       auto d_B_i = B_i.Write();
1196       auto d_X_r = X_r.Read();
1197       auto d_X_i = X_i.Read();
1198       auto d_idx = ess_tdof_list.Read();
1199       MFEM_FORALL(i, n,
1200       {
1201          const int j = d_idx[i];
1202          d_B_r[j] = d_X_r[j];
1203          d_B_i[j] = d_X_i[j];
1204       });
1205       // Modify offdiagonal blocks (imaginary parts of the matrix) to conform
1206       // with standard essential BC treatment
1207       if (A_i.Type() == Operator::Hypre_ParCSR)
1208       {
1209          HypreParMatrix * Ah;
1210          A_i.Get(Ah);
1211          hypre_ParCSRMatrix *Aih = *Ah;
1212 #ifndef HYPRE_USING_CUDA
1213          ess_tdof_list.HostRead();
1214          for (int k = 0; k < n; k++)
1215          {
1216             const int j = ess_tdof_list[k];
1217             Aih->diag->data[Aih->diag->i[j]] = 0.0;
1218          }
1219 #else
1220          Ah->HypreReadWrite();
1221          const int *d_ess_tdof_list =
1222             ess_tdof_list.GetMemory().Read(MemoryClass::DEVICE, n);
1223          const int *d_diag_i = Aih->diag->i;
1224          double *d_diag_data = Aih->diag->data;
1225          CuWrap1D(n, [=] MFEM_DEVICE (int k)
1226          {
1227             const int j = d_ess_tdof_list[k];
1228             d_diag_data[d_diag_i[j]] = 0.0;
1229          });
1230 #endif
1231       }
1232       else
1233       {
1234          A_i.As<ConstrainedOperator>()->SetDiagonalPolicy
1235          (mfem::Operator::DiagonalPolicy::DIAG_ZERO);
1236       }
1237    }
1238 
1239    if (conv == ComplexOperator::BLOCK_SYMMETRIC)
1240    {
1241       B_i *= -1.0;
1242       b_i *= -1.0;
1243    }
1244 
1245    x_r.SyncAliasMemory(x);
1246    x_i.SyncAliasMemory(x);
1247    b_r.SyncAliasMemory(b);
1248    b_i.SyncAliasMemory(b);
1249 
1250    X_r.SyncAliasMemory(X);
1251    X_i.SyncAliasMemory(X);
1252    B_r.SyncAliasMemory(B);
1253    B_i.SyncAliasMemory(B);
1254 
1255    // A = A_r + i A_i
1256    A.Clear();
1257    if ( A_r.Type() == Operator::Hypre_ParCSR ||
1258         A_i.Type() == Operator::Hypre_ParCSR )
1259    {
1260       ComplexHypreParMatrix * A_hyp =
1261          new ComplexHypreParMatrix(A_r.As<HypreParMatrix>(),
1262                                    A_i.As<HypreParMatrix>(),
1263                                    A_r.OwnsOperator(),
1264                                    A_i.OwnsOperator(),
1265                                    conv);
1266       A.Reset<ComplexHypreParMatrix>(A_hyp, true);
1267    }
1268    else
1269    {
1270       ComplexOperator * A_op =
1271          new ComplexOperator(A_r.As<Operator>(),
1272                              A_i.As<Operator>(),
1273                              A_r.OwnsOperator(),
1274                              A_i.OwnsOperator(),
1275                              conv);
1276       A.Reset<ComplexOperator>(A_op, true);
1277    }
1278    A_r.SetOperatorOwner(false);
1279    A_i.SetOperatorOwner(false);
1280 }
1281 
1282 void
FormSystemMatrix(const Array<int> & ess_tdof_list,OperatorHandle & A)1283 ParSesquilinearForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
1284                                       OperatorHandle &A)
1285 {
1286    OperatorHandle A_r, A_i;
1287    if (RealInteg())
1288    {
1289       pblfr->FormSystemMatrix(ess_tdof_list, A_r);
1290    }
1291    if (ImagInteg())
1292    {
1293       pblfi->FormSystemMatrix(ess_tdof_list, A_i);
1294    }
1295    if (!RealInteg() && !ImagInteg())
1296    {
1297       MFEM_ABORT("Both Real and Imaginary part of the Sesquilinear form are empty");
1298    }
1299 
1300    if (RealInteg() && ImagInteg())
1301    {
1302       // Modify offdiagonal blocks (imaginary parts of the matrix) to conform
1303       // with standard essential BC treatment
1304       if ( A_i.Type() == Operator::Hypre_ParCSR )
1305       {
1306          int n = ess_tdof_list.Size();
1307          HypreParMatrix * Ah;
1308          A_i.Get(Ah);
1309          hypre_ParCSRMatrix * Aih = *Ah;
1310          for (int k = 0; k < n; k++)
1311          {
1312             int j = ess_tdof_list[k];
1313             Aih->diag->data[Aih->diag->i[j]] = 0.0;
1314          }
1315       }
1316       else
1317       {
1318          A_i.As<ConstrainedOperator>()->SetDiagonalPolicy
1319          (mfem::Operator::DiagonalPolicy::DIAG_ZERO);
1320       }
1321    }
1322 
1323    // A = A_r + i A_i
1324    A.Clear();
1325    if ( A_r.Type() == Operator::Hypre_ParCSR ||
1326         A_i.Type() == Operator::Hypre_ParCSR )
1327    {
1328       ComplexHypreParMatrix * A_hyp =
1329          new ComplexHypreParMatrix(A_r.As<HypreParMatrix>(),
1330                                    A_i.As<HypreParMatrix>(),
1331                                    A_r.OwnsOperator(),
1332                                    A_i.OwnsOperator(),
1333                                    conv);
1334       A.Reset<ComplexHypreParMatrix>(A_hyp, true);
1335    }
1336    else
1337    {
1338       ComplexOperator * A_op =
1339          new ComplexOperator(A_r.As<Operator>(),
1340                              A_i.As<Operator>(),
1341                              A_r.OwnsOperator(),
1342                              A_i.OwnsOperator(),
1343                              conv);
1344       A.Reset<ComplexOperator>(A_op, true);
1345    }
1346    A_r.SetOperatorOwner(false);
1347    A_i.SetOperatorOwner(false);
1348 }
1349 
1350 void
RecoverFEMSolution(const Vector & X,const Vector & b,Vector & x)1351 ParSesquilinearForm::RecoverFEMSolution(const Vector &X, const Vector &b,
1352                                         Vector &x)
1353 {
1354    ParFiniteElementSpace *pfes = pblfr->ParFESpace();
1355 
1356    const Operator &P = *pfes->GetProlongationMatrix();
1357 
1358    const int vsize  = pfes->GetVSize();
1359    const int tvsize = X.Size() / 2;
1360 
1361    X.Read();
1362    Vector X_r; X_r.MakeRef(const_cast<Vector&>(X), 0, tvsize);
1363    Vector X_i; X_i.MakeRef(const_cast<Vector&>(X), tvsize, tvsize);
1364 
1365    x.Write();
1366    Vector x_r; x_r.MakeRef(x, 0, vsize);
1367    Vector x_i; x_i.MakeRef(x, vsize, vsize);
1368 
1369    // Apply conforming prolongation
1370    P.Mult(X_r, x_r);
1371    P.Mult(X_i, x_i);
1372 
1373    x_r.SyncAliasMemory(x);
1374    x_i.SyncAliasMemory(x);
1375 }
1376 
1377 void
Update(FiniteElementSpace * nfes)1378 ParSesquilinearForm::Update(FiniteElementSpace *nfes)
1379 {
1380    if ( pblfr ) { pblfr->Update(nfes); }
1381    if ( pblfi ) { pblfi->Update(nfes); }
1382 }
1383 
1384 #endif // MFEM_USE_MPI
1385 
1386 }
1387