MFEM  v4.2.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
maxwell_solver.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, 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 "maxwell_solver.hpp"
13 
14 #ifdef MFEM_USE_MPI
15 
16 using namespace std;
17 
18 namespace mfem
19 {
20 
21 using namespace common;
22 
23 namespace electromagnetics
24 {
25 
26 // Used for combining scalar coefficients
27 double prodFunc(double a, double b) { return a * b; }
28 
29 MaxwellSolver::MaxwellSolver(ParMesh & pmesh, int order,
30  double (*eps )(const Vector&),
31  double (*muInv )(const Vector&),
32  double (*sigma )(const Vector&),
33  void (*j_src )(const Vector&, double, Vector&),
34  Array<int> & abcs,
35  Array<int> & dbcs,
36  void (*dEdt_bc )(const Vector&, double, Vector&))
37  : myid_(0),
38  num_procs_(1),
39  order_(order),
40  logging_(1),
41  dtMax_(-1.0),
42  dtScale_(1.0e6),
43  pmesh_(&pmesh),
44  HCurlFESpace_(NULL),
45  HDivFESpace_(NULL),
46  hDivMassMuInv_(NULL),
47  hCurlLosses_(NULL),
48  weakCurlMuInv_(NULL),
49  Curl_(NULL),
50  e_(NULL),
51  b_(NULL),
52  j_(NULL),
53  dedt_(NULL),
54  rhs_(NULL),
55  jd_(NULL),
56  M1Losses_(NULL),
57  M2MuInv_(NULL),
58  NegCurl_(NULL),
59  WeakCurlMuInv_(NULL),
60  E_(NULL),
61  B_(NULL),
62  HD_(NULL),
63  RHS_(NULL),
64  epsCoef_(NULL),
65  muInvCoef_(NULL),
66  sigmaCoef_(NULL),
67  etaInvCoef_(NULL),
68  eCoef_(NULL),
69  bCoef_(NULL),
70  jCoef_(NULL),
71  dEdtBCCoef_(NULL),
72  eps_(eps),
73  muInv_(muInv),
74  sigma_(sigma),
75  j_src_(j_src),
76  dEdt_bc_(dEdt_bc),
77  visit_dc_(NULL)
78 {
79  // Initialize MPI variables
80  MPI_Comm_size(pmesh_->GetComm(), &num_procs_);
81  MPI_Comm_rank(pmesh_->GetComm(), &myid_);
82 
83  // Define compatible parallel finite element spaces on the parallel
84  // mesh. Here we use arbitrary order H1, Nedelec, and Raviart-Thomas finite
85  // elements.
86  HCurlFESpace_ = new ND_ParFESpace(pmesh_,order_,pmesh_->Dimension());
87  HDivFESpace_ = new RT_ParFESpace(pmesh_,order_,pmesh_->Dimension());
88 
89  this->height = HCurlFESpace_->GlobalTrueVSize();
90  this->width = HDivFESpace_->GlobalTrueVSize();
91 
92  // Check for absorbing materials or boundaries
93  lossy_ = abcs.Size() > 0 || sigma_ != NULL;
94 
95  // Require implicit handling of loss terms
96  type = lossy_ ? IMPLICIT : EXPLICIT;
97 
98  // Electric permittivity
99  if ( eps_ == NULL )
100  {
101  epsCoef_ = new ConstantCoefficient(epsilon0_);
102  }
103  else
104  {
105  if ( myid_ == 0 && logging_ > 0 )
106  {
107  cout << "Creating Permittivity Coefficient" << endl;
108  }
109  epsCoef_ = new FunctionCoefficient(eps_);
110  }
111 
112  // Inverse of the magnetic permeability
113  if ( muInv_ == NULL )
114  {
115  muInvCoef_ = new ConstantCoefficient(1.0/mu0_);
116  }
117  else
118  {
119  if ( myid_ == 0 && logging_ > 0 )
120  {
121  cout << "Creating Permeability Coefficient" << endl;
122  }
123  muInvCoef_ = new FunctionCoefficient(muInv_);
124  }
125 
126  // Electric conductivity
127  if ( sigma_ != NULL )
128  {
129  if ( myid_ == 0 && logging_ > 0 )
130  {
131  cout << "Creating Conductivity Coefficient" << endl;
132  }
133  sigmaCoef_ = new FunctionCoefficient(sigma_);
134  }
135 
136  // Impedance of free space
137  if ( abcs.Size() > 0 )
138  {
139  if ( myid_ == 0 && logging_ > 0 )
140  {
141  cout << "Creating Admittance Coefficient" << endl;
142  }
143 
144  AttrToMarker(pmesh.bdr_attributes.Max(), abcs, abc_marker_);
145  etaInvCoef_ = new ConstantCoefficient(sqrt(epsilon0_/mu0_));
146  }
147 
148  // Electric Field Boundary Condition
149  if ( dbcs.Size() > 0 )
150  {
151  if ( myid_ == 0 && logging_ > 0 )
152  {
153  cout << "Configuring Dirichlet BC" << endl;
154  }
155 
156  AttrToMarker(pmesh.bdr_attributes.Max(), dbcs, dbc_marker_);
157  HCurlFESpace_->GetEssentialTrueDofs(dbc_marker_, dbc_dofs_);
158 
159  if ( dEdt_bc_ != NULL )
160  {
161  dEdtBCCoef_ = new VectorFunctionCoefficient(3,dEdt_bc_);
162  }
163  else
164  {
165  Vector ebc(3); ebc = 0.0;
166  dEdtBCCoef_ = new VectorConstantCoefficient(ebc);
167  }
168  }
169 
170  // Bilinear Forms
171  if ( myid_ == 0 && logging_ > 0 )
172  {
173  cout << "Creating H(Div) Mass Operator" << endl;
174  }
175  hDivMassMuInv_ = new ParBilinearForm(HDivFESpace_);
176  hDivMassMuInv_->AddDomainIntegrator(new VectorFEMassIntegrator(*muInvCoef_));
177 
178  if ( myid_ == 0 && logging_ > 0 )
179  {
180  cout << "Creating Weak Curl Operator" << endl;
181  }
182  weakCurlMuInv_ = new ParMixedBilinearForm(HDivFESpace_,HCurlFESpace_);
183  weakCurlMuInv_->AddDomainIntegrator(
184  new MixedVectorWeakCurlIntegrator(*muInvCoef_));
185 
186  // Assemble Matrices
187  hDivMassMuInv_->Assemble();
188  weakCurlMuInv_->Assemble();
189 
190  hDivMassMuInv_->Finalize();
191  weakCurlMuInv_->Finalize();
192 
193  if ( sigmaCoef_ || etaInvCoef_ )
194  {
195  if ( myid_ == 0 && logging_ > 0 )
196  {
197  cout << "Creating H(Curl) Loss Operator" << endl;
198  }
199  hCurlLosses_ = new ParBilinearForm(HCurlFESpace_);
200  if ( sigmaCoef_ )
201  {
202  if ( myid_ == 0 && logging_ > 0 )
203  {
204  cout << "Adding domain integrator for conductive regions" << endl;
205  }
206  hCurlLosses_->AddDomainIntegrator(
207  new VectorFEMassIntegrator(*sigmaCoef_));
208  }
209  if ( etaInvCoef_ )
210  {
211  if ( myid_ == 0 && logging_ > 0 )
212  {
213  cout << "Adding boundary integrator for absorbing boundary" << endl;
214  }
215  hCurlLosses_->AddBoundaryIntegrator(
216  new VectorFEMassIntegrator(*etaInvCoef_), abc_marker_);
217  }
218  hCurlLosses_->Assemble();
219  hCurlLosses_->Finalize();
220  M1Losses_ = hCurlLosses_->ParallelAssemble();
221  }
222 
223  // Create Linear Algebra Matrices
224  M2MuInv_ = hDivMassMuInv_->ParallelAssemble();
225  WeakCurlMuInv_ = weakCurlMuInv_->ParallelAssemble();
226 
227  if ( myid_ == 0 && logging_ > 0 )
228  {
229  cout << "Creating discrete curl operator" << endl;
230  }
231  Curl_ = new ParDiscreteCurlOperator(HCurlFESpace_, HDivFESpace_);
232  Curl_->Assemble();
233  Curl_->Finalize();
234  NegCurl_ = Curl_->ParallelAssemble();
235  // Beware this modifies the matrix stored within the Curl_ object.
236  *NegCurl_ *= -1.0;
237 
238  // Build grid functions
239  e_ = new ParGridFunction(HCurlFESpace_);
240  dedt_ = new ParGridFunction(HCurlFESpace_);
241  rhs_ = new ParGridFunction(HCurlFESpace_);
242  b_ = new ParGridFunction(HDivFESpace_);
243 
244  E_ = e_->ParallelProject();
245  B_ = b_->ParallelProject();
246 
247  HD_ = new HypreParVector(HDivFESpace_);
248  RHS_ = new HypreParVector(HCurlFESpace_);
249 
250  // Initialize dedt to zero
251  *dedt_ = 0.0;
252 
253  if ( j_src_)
254  {
255  if ( myid_ == 0 && logging_ > 0 )
256  {
257  cout << "Creating Current Source" << endl;
258  }
259  jCoef_ = new VectorFunctionCoefficient(3,j_src_);
260 
261  j_ = new ParGridFunction(HCurlFESpace_);
262  j_->ProjectCoefficient(*jCoef_);
263 
264  jd_ = new ParLinearForm(HCurlFESpace_);
266  jd_->Assemble();
267  }
268 
269  dtMax_ = GetMaximumTimeStep();
270 }
271 
273 {
274  delete epsCoef_;
275  delete muInvCoef_;
276  delete etaInvCoef_;
277  delete jCoef_;
278  delete dEdtBCCoef_;
279 
280  delete E_;
281  delete B_;
282  delete HD_;
283  delete RHS_;
284 
285  delete e_;
286  delete b_;
287  delete j_;
288  delete dedt_;
289  delete rhs_;
290  delete jd_;
291 
292  delete Curl_;
293 
294  delete M1Losses_;
295  delete M2MuInv_;
296  delete NegCurl_;
297  delete WeakCurlMuInv_;
298 
299  delete hDivMassMuInv_;
300  delete hCurlLosses_;
301  delete weakCurlMuInv_;
302 
303  delete HCurlFESpace_;
304  delete HDivFESpace_;
305 
306  map<int, ParBilinearForm*>::iterator mit1;
307  for (mit1=a1_.begin(); mit1!=a1_.end(); mit1++)
308  {
309  int i = mit1->first;
310  delete pcg_[i];
311  delete diagScale_[i];
312  delete A1_[i];
313  delete a1_[i];
314  }
315 
316  map<int, Coefficient*>::iterator mit2;
317  for (mit2=dtCoef_.begin(); mit2!=dtCoef_.end(); mit2++)
318  {
319  delete mit2->second;
320  }
321  for (mit2=dtSigmaCoef_.begin(); mit2!=dtSigmaCoef_.end(); mit2++)
322  {
323  delete mit2->second;
324  }
325  for (mit2=dtEtaInvCoef_.begin(); mit2!=dtEtaInvCoef_.end(); mit2++)
326  {
327  delete mit2->second;
328  }
329 
330  map<string, socketstream*>::iterator mit3;
331  for (mit3=socks_.begin(); mit3!=socks_.end(); mit3++)
332  {
333  delete mit3->second;
334  }
335 }
336 
337 HYPRE_Int
339 {
340  return HCurlFESpace_->GlobalTrueVSize();
341 }
342 
343 void
345 {
346  HYPRE_Int size_nd = HCurlFESpace_->GlobalTrueVSize();
347  HYPRE_Int size_rt = HDivFESpace_->GlobalTrueVSize();
348  if ( myid_ == 0 )
349  {
350  cout << "Number of H(Curl) unknowns: " << size_nd << endl;
351  cout << "Number of H(Div) unknowns: " << size_rt << endl << flush;
352  }
353 }
354 
355 void
357 {
358  eCoef_ = &EFieldCoef;
359  e_->ProjectCoefficient(EFieldCoef);
360  e_->ParallelProject(*E_);
361 }
362 
363 void
365 {
366  bCoef_ = &BFieldCoef;
367  b_->ProjectCoefficient(BFieldCoef);
368  b_->ParallelProject(*B_);
369 }
370 
371 void
372 MaxwellSolver::Mult(const Vector &B, Vector &dEdt) const
373 {
374  implicitSolve(0.0, B, dEdt);
375 }
376 
377 void
378 MaxwellSolver::ImplicitSolve(double dt, const Vector &B, Vector &dEdt)
379 {
380  implicitSolve(dt, B, dEdt);
381 }
382 
383 void
384 MaxwellSolver::setupSolver(const int idt, const double dt) const
385 {
386  if ( pcg_.find(idt) == pcg_.end() )
387  {
388  if ( myid_ == 0 && logging_ > 0 )
389  {
390  cout << "Creating implicit operator for dt = " << dt << endl;
391  }
392 
393  a1_[idt] = new ParBilinearForm(HCurlFESpace_);
394  a1_[idt]->AddDomainIntegrator(
395  new VectorFEMassIntegrator(epsCoef_));
396 
397  if ( idt != 0 )
398  {
399  dtCoef_[idt] = new ConstantCoefficient(0.5 * dt);
400  if ( sigmaCoef_ )
401  {
402  dtSigmaCoef_[idt] = new TransformedCoefficient(dtCoef_[idt],
403  sigmaCoef_,
404  prodFunc);
405  a1_[idt]->AddDomainIntegrator(
406  new VectorFEMassIntegrator(dtSigmaCoef_[idt]));
407  }
408  if ( etaInvCoef_ )
409  {
410  dtEtaInvCoef_[idt] = new TransformedCoefficient(dtCoef_[idt],
411  etaInvCoef_,
412  prodFunc);
413  a1_[idt]->AddBoundaryIntegrator(
414  new VectorFEMassIntegrator(dtEtaInvCoef_[idt]),
415  const_cast<Array<int>&>(abc_marker_));
416  }
417  }
418 
419  a1_[idt]->Assemble();
420  a1_[idt]->Finalize();
421  A1_[idt] = a1_[idt]->ParallelAssemble();
422 
423  diagScale_[idt] = new HypreDiagScale(*A1_[idt]);
424  pcg_[idt] = new HyprePCG(*A1_[idt]);
425  pcg_[idt]->SetTol(1.0e-12);
426  pcg_[idt]->SetMaxIter(200);
427  pcg_[idt]->SetPrintLevel(0);
428  pcg_[idt]->SetPreconditioner(*diagScale_[idt]);
429  }
430 }
431 
432 void
433 MaxwellSolver::implicitSolve(double dt, const Vector &B, Vector &dEdt) const
434 {
435  int idt = hCurlLosses_ ? ((int)(dtScale_ * dt / dtMax_)) : 0;
436 
437  b_->Distribute(B);
438  weakCurlMuInv_->Mult(*b_, *rhs_);
439 
440  if ( hCurlLosses_ )
441  {
442  e_->Distribute(*E_);
443  hCurlLosses_->AddMult(*e_, *rhs_, -1.0);
444  }
445 
446  if ( jd_ )
447  {
448  jCoef_->SetTime(t); // 't' is member data from mfem::TimeDependentOperator
449  jd_->Assemble();
450  *rhs_ -= *jd_;
451  }
452 
453  if ( dEdtBCCoef_ )
454  {
455  dEdtBCCoef_->SetTime(t);
456  dedt_->ProjectBdrCoefficientTangent(*dEdtBCCoef_,
457  const_cast<Array<int>&>(dbc_marker_));
458  }
459 
460  // Create objects and matrices for solving with the given time step
461  setupSolver(idt, dt);
462 
463  // Apply essential BCs and determine true DoFs for the right hand side
464  a1_[idt]->FormLinearSystem(dbc_dofs_, *dedt_, *rhs_, *A1_[idt], dEdt, *RHS_);
465 
466  // Solve for the time derivative of the electric field (true DoFs)
467  pcg_[idt]->Mult(*RHS_, dEdt);
468 
469  // Distribute shared DoFs to relevant processors
470  a1_[idt]->RecoverFEMSolution(dEdt, *rhs_, *dedt_);
471 }
472 
473 void
475 {
476  e_->Distribute(*E_);
477  b_->Distribute(*B_);
478 }
479 
480 double
482 {
483  if ( dtMax_ > 0.0 )
484  {
485  return dtMax_;
486  }
487 
488  HypreParVector * v0 = new HypreParVector(HCurlFESpace_);
489  HypreParVector * v1 = new HypreParVector(HCurlFESpace_);
490  HypreParVector * u0 = new HypreParVector(HDivFESpace_);
491 
492  v0->Randomize(1234);
493 
494  int iter = 0, nstep = 20;
495  double dt0 = 1.0, dt1 = 1.0, change = 1.0, ptol = 0.001;
496 
497  // Create Solver assuming no loss operators
498  setupSolver(0, 0.0);
499 
500  // Use power method to approximate the largest eigenvalue of the update
501  // operator.
502  while ( iter < nstep && change > ptol )
503  {
504  double normV0 = InnerProduct(*v0,*v0);
505  *v0 /= sqrt(normV0);
506 
507  NegCurl_->Mult(*v0,*u0);
508  M2MuInv_->Mult(*u0,*HD_);
509  NegCurl_->MultTranspose(*HD_,*RHS_);
510 
511  pcg_[0]->Mult(*RHS_,*v1);
512 
513  double lambda = InnerProduct(*v0,*v1);
514  dt1 = 2.0/sqrt(lambda);
515  change = fabs((dt1-dt0)/dt0);
516  dt0 = dt1;
517 
518  if ( myid_ == 0 && logging_ > 1 )
519  {
520  cout << iter << ": " << dt0 << " " << change << endl;
521  }
522 
523  std::swap(v0, v1);
524 
525  iter++;
526  }
527 
528  delete v0;
529  delete v1;
530  delete u0;
531 
532  return dt0;
533 }
534 
535 double
537 {
538  double energy = 0.0;
539 
540  A1_[0]->Mult(*E_,*RHS_);
541  M2MuInv_->Mult(*B_,*HD_);
542 
543  energy = InnerProduct(*E_,*RHS_) + InnerProduct(*B_,*HD_);
544 
545  return 0.5 * energy;
546 }
547 
548 void
550 {
551  visit_dc_ = &visit_dc;
552 
553  visit_dc.RegisterField("E", e_);
554  visit_dc.RegisterField("B", b_);
555  if ( j_ )
556  {
557  visit_dc.RegisterField("J", j_);
558  }
559 }
560 
561 void
563 {
564  if ( visit_dc_ )
565  {
566  if ( myid_ == 0 && logging_ > 1 )
567  { cout << "Writing VisIt files ..." << flush; }
568 
569  if ( j_ )
570  {
571  jCoef_->SetTime(t);
572  j_->ProjectCoefficient(*jCoef_);
573  }
574 
575  visit_dc_->SetCycle(it);
576  visit_dc_->SetTime(t);
577  visit_dc_->Save();
578 
579  if ( myid_ == 0 && logging_ > 1 ) { cout << " " << endl << flush; }
580  }
581 }
582 
583 void
585 {
586  if ( myid_ == 0 && logging_ > 0 )
587  { cout << "Opening GLVis sockets." << endl << flush; }
588 
589  socks_["E"] = new socketstream;
590  socks_["E"]->precision(8);
591 
592  socks_["B"] = new socketstream;
593  socks_["B"]->precision(8);
594 
595  if ( j_ )
596  {
597  socks_["J"] = new socketstream;
598  socks_["J"]->precision(8);
599  }
600 
601  if ( myid_ == 0 && logging_ > 0 )
602  { cout << "GLVis sockets open." << endl << flush; }
603 }
604 
605 void
607 {
608  if ( myid_ == 0 && logging_ > 1 )
609  { cout << "Sending data to GLVis ..." << flush; }
610 
611  char vishost[] = "localhost";
612  int visport = 19916;
613 
614  int Wx = 0, Wy = 0; // window position
615  int Ww = 350, Wh = 350; // window size
616  int offx = Ww+10, offy = Wh+45; // window offsets
617 
618  VisualizeField(*socks_["E"], vishost, visport,
619  *e_, "Electric Field (E)", Wx, Wy, Ww, Wh);
620  Wx += offx;
621 
622  VisualizeField(*socks_["B"], vishost, visport,
623  *b_, "Magnetic Flux Density (B)", Wx, Wy, Ww, Wh);
624 
625  if ( j_ )
626  {
627  Wx = 0;
628  Wy += offy;
629 
630  jCoef_->SetTime(t); // Is member data from mfem::TimeDependentOperator
631  j_->ProjectCoefficient(*jCoef_);
632 
633  VisualizeField(*socks_["J"], vishost, visport,
634  *j_, "Current Density (J)", Wx, Wy, Ww, Wh);
635  }
636  if ( myid_ == 0 && logging_ > 1 ) { cout << " " << flush; }
637 }
638 
639 } // namespace electromagnetics
640 
641 } // namespace mfem
642 
643 #endif // MFEM_USE_MPI
A coefficient that depends on 1 or 2 parent coefficients and a transformation rule represented by a C...
int Size() const
Return the logical size of the array.
Definition: array.hpp:124
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
Base class for vector Coefficients that optionally depend on time and space.
A coefficient that is constant across space and time.
Definition: coefficient.hpp:78
void AttrToMarker(int max_attr, const Array< int > &attrs, Array< int > &marker)
Convert a set of attribute numbers to a marker array.
Type type
Describes the form of the TimeDependentOperator.
Definition: operator.hpp:292
HYPRE_Int MultTranspose(HypreParVector &x, HypreParVector &y, double alpha=1.0, double beta=0.0)
Computes y = alpha * A^t * x + beta * y.
Definition: hypre.cpp:1045
virtual void AddMult(const Vector &x, Vector &y, const double a=1.0) const
Add the matrix vector multiple to a vector: .
double prodFunc(double a, double b)
Vector coefficient that is constant in space and time.
int offx
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:79
int Wx
virtual void Save()
Save the collection and a VisIt root file.
virtual void ProjectCoefficient(Coefficient &coeff)
Definition: pgridfunc.cpp:474
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
HYPRE_Int Mult(HypreParVector &x, HypreParVector &y, double alpha=1.0, double beta=0.0)
Computes y = alpha * A * x + beta * y.
Definition: hypre.cpp:969
void SetInitialBField(VectorCoefficient &BFieldCoef)
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
double t
Current time.
Definition: operator.hpp:291
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P_test^t A P_trial.
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds a domain integrator. Assumes ownership of bfi.
void ImplicitSolve(const double dt, const Vector &x, Vector &k)
Solve the equation: k = f(x + dt k, t), for the unknown k at the current time t.
Class for parallel linear form.
Definition: plinearform.hpp:26
constexpr char vishost[]
void RegisterVisItFields(VisItDataCollection &visit_dc)
double b
Definition: lissajous.cpp:42
This type assumes F(x,k,t) = k, i.e. k = f(x,t) = G(x,t).
Definition: operator.hpp:272
constexpr int visport
Data collection with VisIt I/O routines.
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
void Assemble(int skip_zeros=1)
double muInv(const Vector &x)
Definition: maxwell.cpp:96
HYPRE_Int GlobalTrueVSize() const
Definition: pfespace.hpp:251
void Assemble(int skip_zeros=1)
Assemble the local matrix.
HYPRE_Int Randomize(HYPRE_Int seed)
Set random values.
Definition: hypre.cpp:166
This is the most general type, no assumptions on F and G.
Definition: operator.hpp:273
int Dimension() const
Definition: mesh.hpp:788
void SetTime(double t)
Set physical time (for time-dependent simulations)
A general vector function coefficient.
HypreParMatrix * ParallelAssemble() const
Returns the matrix &quot;assembled&quot; on the true dofs.
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:70
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:39
void Mult(const Vector &B, Vector &dEdt) const
Perform the action of the operator: y = k = f(x, t), where k solves the algebraic equation F(x...
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:201
void j_src(const Vector &x, double t, Vector &j)
Definition: maxwell.cpp:108
MPI_Comm GetComm() const
Definition: pmesh.hpp:236
void SetTime(double t)
Set the time for time dependent coefficients.
void Distribute(const Vector *tv)
Definition: pgridfunc.cpp:132
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection and update the root file.
double a
Definition: lissajous.cpp:41
double InnerProduct(HypreParVector *x, HypreParVector *y)
Definition: hypre.cpp:194
Class for parallel bilinear form using different test and trial FE spaces.
int Wh
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
void AddBoundaryIntegrator(BilinearFormIntegrator *bfi)
Adds new Boundary Integrator. Assumes ownership of bfi.
virtual void ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff, Array< int > &bdr_attr)
Project the tangential components of the given VectorCoefficient on the boundary. Only boundary attri...
Definition: pgridfunc.cpp:620
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
Class for parallel bilinear form.
for VectorFiniteElements (Nedelec, Raviart-Thomas)
Definition: lininteg.hpp:260
A general function coefficient.
Vector data type.
Definition: vector.hpp:51
void VisualizeField(socketstream &sock, const char *vishost, int visport, GridFunction &gf, const char *title, int x, int y, int w, int h, const char *keys, bool vec)
Definition: fem_extras.cpp:94
int offy
int Wy
int Ww
Class for parallel grid function.
Definition: pgridfunc.hpp:32
virtual void Assemble(int skip_zeros=1)
Construct the internal matrix representation of the discrete linear operator.
Class for parallel meshes.
Definition: pmesh.hpp:32
virtual void Mult(const Vector &x, Vector &y) const
Matrix multiplication: .
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28
void SetInitialEField(VectorCoefficient &EFieldCoef)
void ParallelProject(Vector &tv) const
Returns the vector restricted to the true dofs.
Definition: pgridfunc.cpp:171
double sigma(const Vector &x)
Definition: maxwell.cpp:102