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