MFEM  v4.6.0
Finite element discretization library
tesla_solver.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 "tesla_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 TeslaSolver::TeslaSolver(ParMesh & pmesh, int order,
27  Array<int> & kbcs,
28  Array<int> & vbcs, Vector & vbcv,
29  Coefficient & muInvCoef,
30  void (*a_bc )(const Vector&, Vector&),
31  void (*j_src)(const Vector&, Vector&),
32  void (*m_src)(const Vector&, Vector&))
33  : myid_(0),
34  num_procs_(1),
35  order_(order),
36  pmesh_(&pmesh),
37  visit_dc_(NULL),
38  H1FESpace_(NULL),
39  HCurlFESpace_(NULL),
40  HDivFESpace_(NULL),
41  curlMuInvCurl_(NULL),
42  hCurlMass_(NULL),
43  hDivHCurlMuInv_(NULL),
44  weakCurlMuInv_(NULL),
45  grad_(NULL),
46  curl_(NULL),
47  a_(NULL),
48  b_(NULL),
49  h_(NULL),
50  jr_(NULL),
51  j_(NULL),
52  k_(NULL),
53  m_(NULL),
54  bd_(NULL),
55  jd_(NULL),
56  DivFreeProj_(NULL),
57  SurfCur_(NULL),
58  muInvCoef_(&muInvCoef),
59  aBCCoef_(NULL),
60  jCoef_(NULL),
61  mCoef_(NULL),
62  a_bc_(a_bc),
63  j_src_(j_src),
64  m_src_(m_src)
65 {
66  // Initialize MPI variables
67  MPI_Comm_size(pmesh_->GetComm(), &num_procs_);
68  MPI_Comm_rank(pmesh_->GetComm(), &myid_);
69 
70  // Define compatible parallel finite element spaces on the parallel
71  // mesh. Here we use arbitrary order H1, Nedelec, and Raviart-Thomas finite
72  // elements.
73  H1FESpace_ = new H1_ParFESpace(pmesh_,order,pmesh_->Dimension());
74  HCurlFESpace_ = new ND_ParFESpace(pmesh_,order,pmesh_->Dimension());
75  HDivFESpace_ = new RT_ParFESpace(pmesh_,order,pmesh_->Dimension());
76 
77  int irOrder = H1FESpace_->GetElementTransformation(0)->OrderW()
78  + 2 * order;
79  int geom = H1FESpace_->GetFE(0)->GetGeomType();
80  const IntegrationRule * ir = &IntRules.Get(geom, irOrder);
81 
82  // Select surface attributes for Dirichlet BCs
83  ess_bdr_.SetSize(pmesh.bdr_attributes.Max());
84  non_k_bdr_.SetSize(pmesh.bdr_attributes.Max());
85  ess_bdr_ = 1; // All outer surfaces
86  non_k_bdr_ = 1; // Surfaces without applied surface currents
87 
88  for (int i=0; i<kbcs.Size(); i++)
89  {
90  non_k_bdr_[kbcs[i]-1] = 0;
91  }
92 
93  // Setup various coefficients
94 
95  // Vector Potential on the outer surface
96  if ( a_bc_ == NULL )
97  {
98  Vector Zero(3);
99  Zero = 0.0;
100  aBCCoef_ = new VectorConstantCoefficient(Zero);
101  }
102  else
103  {
104  aBCCoef_ = new VectorFunctionCoefficient(pmesh_->SpaceDimension(),
105  *a_bc_);
106  }
107 
108  // Volume Current Density
109  if ( j_src_ != NULL )
110  {
111  jCoef_ = new VectorFunctionCoefficient(pmesh_->SpaceDimension(),
112  j_src_);
113  }
114 
115  // Magnetization
116  if ( m_src_ != NULL )
117  {
118  mCoef_ = new VectorFunctionCoefficient(pmesh_->SpaceDimension(),
119  m_src_);
120  }
121 
122  // Bilinear Forms
123  curlMuInvCurl_ = new ParBilinearForm(HCurlFESpace_);
124  curlMuInvCurl_->AddDomainIntegrator(new CurlCurlIntegrator(*muInvCoef_));
125 
126  BilinearFormIntegrator * hCurlMassInteg = new VectorFEMassIntegrator;
127  hCurlMassInteg->SetIntRule(ir);
128  hCurlMass_ = new ParBilinearForm(HCurlFESpace_);
129  hCurlMass_->AddDomainIntegrator(hCurlMassInteg);
130 
131  BilinearFormIntegrator * hDivHCurlInteg =
132  new VectorFEMassIntegrator(*muInvCoef_);
133  hDivHCurlInteg->SetIntRule(ir);
134  hDivHCurlMuInv_ = new ParMixedBilinearForm(HDivFESpace_, HCurlFESpace_);
135  hDivHCurlMuInv_->AddDomainIntegrator(hDivHCurlInteg);
136 
137  // Discrete Curl operator
138  curl_ = new ParDiscreteCurlOperator(HCurlFESpace_, HDivFESpace_);
139 
140  // Build grid functions
141  a_ = new ParGridFunction(HCurlFESpace_);
142  b_ = new ParGridFunction(HDivFESpace_);
143  h_ = new ParGridFunction(HCurlFESpace_);
144  bd_ = new ParGridFunction(HCurlFESpace_);
145  jd_ = new ParGridFunction(HCurlFESpace_);
146 
147  if ( jCoef_ || kbcs.Size() > 0 )
148  {
149  grad_ = new ParDiscreteGradOperator(H1FESpace_, HCurlFESpace_);
150  }
151  if ( jCoef_ )
152  {
153  jr_ = new ParGridFunction(HCurlFESpace_);
154  j_ = new ParGridFunction(HCurlFESpace_);
155  DivFreeProj_ = new DivergenceFreeProjector(*H1FESpace_, *HCurlFESpace_,
156  irOrder, NULL, NULL, grad_);
157  }
158 
159  if ( kbcs.Size() > 0 )
160  {
161  k_ = new ParGridFunction(HCurlFESpace_);
162 
163  // Object to solve the subproblem of computing surface currents
164  SurfCur_ = new SurfaceCurrent(*H1FESpace_, *grad_,
165  kbcs, vbcs, vbcv);
166  }
167 
168  if ( mCoef_ )
169  {
170  m_ = new ParGridFunction(HDivFESpace_);
171 
172  weakCurlMuInv_ = new ParMixedBilinearForm(HDivFESpace_, HCurlFESpace_);
173  weakCurlMuInv_->AddDomainIntegrator(
174  new VectorFECurlIntegrator(*muInvCoef_));
175  }
176 }
177 
179 {
180  delete jCoef_;
181  delete mCoef_;
182  delete aBCCoef_;
183 
184  delete DivFreeProj_;
185  delete SurfCur_;
186 
187  delete a_;
188  delete b_;
189  delete h_;
190  delete jr_;
191  delete j_;
192  delete k_;
193  delete m_;
194  delete bd_;
195  delete jd_;
196 
197  delete grad_;
198  delete curl_;
199 
200  delete curlMuInvCurl_;
201  delete hCurlMass_;
202  delete hDivHCurlMuInv_;
203  delete weakCurlMuInv_;
204 
205  delete H1FESpace_;
206  delete HCurlFESpace_;
207  delete HDivFESpace_;
208 
209  map<string,socketstream*>::iterator mit;
210  for (mit=socks_.begin(); mit!=socks_.end(); mit++)
211  {
212  delete mit->second;
213  }
214 }
215 
218 {
219  return HCurlFESpace_->GlobalTrueVSize();
220 }
221 
222 void
224 {
225  HYPRE_BigInt size_h1 = H1FESpace_->GlobalTrueVSize();
226  HYPRE_BigInt size_nd = HCurlFESpace_->GlobalTrueVSize();
227  HYPRE_BigInt size_rt = HDivFESpace_->GlobalTrueVSize();
228  if (myid_ == 0)
229  {
230  cout << "Number of H1 unknowns: " << size_h1 << endl;
231  cout << "Number of H(Curl) unknowns: " << size_nd << endl;
232  cout << "Number of H(Div) unknowns: " << size_rt << endl;
233  }
234 }
235 
236 void
238 {
239  if (myid_ == 0) { cout << "Assembling ..." << flush; }
240 
241  curlMuInvCurl_->Assemble();
242  curlMuInvCurl_->Finalize();
243 
244  hDivHCurlMuInv_->Assemble();
245  hDivHCurlMuInv_->Finalize();
246 
247  hCurlMass_->Assemble();
248  hCurlMass_->Finalize();
249 
250  curl_->Assemble();
251  curl_->Finalize();
252 
253  if ( grad_ )
254  {
255  grad_->Assemble();
256  grad_->Finalize();
257  }
258  if ( weakCurlMuInv_ )
259  {
260  weakCurlMuInv_->Assemble();
261  weakCurlMuInv_->Finalize();
262  }
263 
264  if (myid_ == 0) { cout << " done." << endl; }
265 }
266 
267 void
269 {
270  if (myid_ == 0) { cout << "Updating ..." << endl; }
271 
272  // Inform the spaces that the mesh has changed
273  // Note: we don't need to interpolate any GridFunctions on the new mesh
274  // so we pass 'false' to skip creation of any transformation matrices.
275  H1FESpace_->Update(false);
276  HCurlFESpace_->Update(false);
277  HDivFESpace_->Update(false);
278 
279  HCurlFESpace_->GetEssentialTrueDofs(ess_bdr_, ess_bdr_tdofs_);
280 
281  // Inform the grid functions that the space has changed.
282  a_->Update();
283  h_->Update();
284  b_->Update();
285  bd_->Update();
286  jd_->Update();
287  if ( jr_ ) { jr_->Update(); }
288  if ( j_ ) { j_->Update(); }
289  if ( k_ ) { k_->Update(); }
290  if ( m_ ) { m_->Update(); }
291 
292  // Inform the bilinear forms that the space has changed.
293  curlMuInvCurl_->Update();
294  hCurlMass_->Update();
295  hDivHCurlMuInv_->Update();
296  if ( weakCurlMuInv_ ) { weakCurlMuInv_->Update(); }
297 
298  // Inform the other objects that the space has changed.
299  curl_->Update();
300  if ( grad_ ) { grad_->Update(); }
301  if ( DivFreeProj_ ) { DivFreeProj_->Update(); }
302  if ( SurfCur_ ) { SurfCur_->Update(); }
303 }
304 
305 void
307 {
308  if (myid_ == 0) { cout << "Running solver ... " << endl; }
309 
310  // Initialize the magnetic vector potential with its boundary conditions
311  *a_ = 0.0;
312 
313  // Apply surface currents if available
314  if ( k_ )
315  {
316  SurfCur_->ComputeSurfaceCurrent(*k_);
317  *a_ = *k_;
318  }
319 
320  // Apply uniform B boundary condition on remaining surfaces
321  a_->ProjectBdrCoefficientTangent(*aBCCoef_, non_k_bdr_);
322 
323  // Initialize the RHS vector to zero
324  *jd_ = 0.0;
325 
326  // Initialize the volumetric current density
327  if ( jr_ )
328  {
329  jr_->ProjectCoefficient(*jCoef_);
330 
331  // Compute the discretely divergence-free portion of jr_
332  DivFreeProj_->Mult(*jr_, *j_);
333 
334  // Compute the dual of j_
335  hCurlMass_->AddMult(*j_, *jd_);
336  }
337 
338  // Initialize the Magnetization
339  if ( m_ )
340  {
341  m_->ProjectCoefficient(*mCoef_);
342  weakCurlMuInv_->AddMult(*m_, *jd_, mu0_);
343  }
344 
345  // Apply Dirichlet BCs to matrix and right hand side and otherwise
346  // prepare the linear system
347  HypreParMatrix CurlMuInvCurl;
348  HypreParVector A(HCurlFESpace_);
349  HypreParVector RHS(HCurlFESpace_);
350 
351  curlMuInvCurl_->FormLinearSystem(ess_bdr_tdofs_, *a_, *jd_, CurlMuInvCurl,
352  A, RHS);
353 
354  // Define and apply a parallel PCG solver for AX=B with the AMS
355  // preconditioner from hypre.
356  HypreAMS ams(CurlMuInvCurl, HCurlFESpace_);
357  ams.SetSingularProblem();
358 
359  HyprePCG pcg (CurlMuInvCurl);
360  pcg.SetTol(1e-12);
361  pcg.SetMaxIter(50);
362  pcg.SetPrintLevel(2);
363  pcg.SetPreconditioner(ams);
364  pcg.Mult(RHS, A);
365 
366  // Extract the parallel grid function corresponding to the finite
367  // element approximation A. This is the local solution on each
368  // processor.
369  curlMuInvCurl_->RecoverFEMSolution(A, *jd_, *a_);
370 
371  // Compute the negative Gradient of the solution vector. This is
372  // the magnetic field corresponding to the scalar potential
373  // represented by phi.
374  curl_->Mult(*a_, *b_);
375 
376  // Compute magnetic field (H) from B and M
377  if (myid_ == 0) { cout << "Computing H ... " << flush; }
378 
379  hDivHCurlMuInv_->Mult(*b_, *bd_);
380  if ( m_ )
381  {
382  hDivHCurlMuInv_->AddMult(*m_, *bd_, -1.0 * mu0_);
383  }
384 
385  HypreParMatrix MassHCurl;
386  Vector BD, H;
387 
388  Array<int> dbc_dofs_h;
389  hCurlMass_->FormLinearSystem(dbc_dofs_h, *h_, *bd_, MassHCurl, H, BD);
390 
391  HyprePCG pcgM(MassHCurl);
392  pcgM.SetTol(1e-12);
393  pcgM.SetMaxIter(500);
394  pcgM.SetPrintLevel(0);
395  HypreDiagScale diagM;
396  pcgM.SetPreconditioner(diagM);
397  pcgM.Mult(BD, H);
398 
399  hCurlMass_->RecoverFEMSolution(H, *bd_, *h_);
400 
401  if (myid_ == 0) { cout << "done." << flush; }
402 
403  if (myid_ == 0) { cout << " Solver done. " << endl; }
404 }
405 
406 void
408 {
409  if (myid_ == 0) { cout << "Estimating Error ... " << flush; }
410 
411  // Space for the discontinuous (original) flux
412  CurlCurlIntegrator flux_integrator(*muInvCoef_);
413  RT_FECollection flux_fec(order_-1, pmesh_->SpaceDimension());
414  ParFiniteElementSpace flux_fes(pmesh_, &flux_fec);
415 
416  // Space for the smoothed (conforming) flux
417  int norm_p = 1;
418  ND_FECollection smooth_flux_fec(order_, pmesh_->Dimension());
419  ParFiniteElementSpace smooth_flux_fes(pmesh_, &smooth_flux_fec);
420 
421  L2ZZErrorEstimator(flux_integrator, *a_,
422  smooth_flux_fes, flux_fes, errors, norm_p);
423 
424  if (myid_ == 0) { cout << "done." << endl; }
425 }
426 
427 void
429 {
430  visit_dc_ = &visit_dc;
431 
432  visit_dc.RegisterField("A", a_);
433  visit_dc.RegisterField("B", b_);
434  visit_dc.RegisterField("H", h_);
435  if ( j_ ) { visit_dc.RegisterField("J", j_); }
436  if ( k_ ) { visit_dc.RegisterField("K", k_); }
437  if ( m_ ) { visit_dc.RegisterField("M", m_); }
438  if ( SurfCur_ ) { visit_dc.RegisterField("Psi", SurfCur_->GetPsi()); }
439 }
440 
441 void
443 {
444  if ( visit_dc_ )
445  {
446  if (myid_ == 0) { cout << "Writing VisIt files ..." << flush; }
447 
448  HYPRE_BigInt prob_size = this->GetProblemSize();
449  visit_dc_->SetCycle(it);
450  visit_dc_->SetTime(prob_size);
451  visit_dc_->Save();
452 
453  if (myid_ == 0) { cout << " done." << endl; }
454  }
455 }
456 
457 void
459 {
460  if ( myid_ == 0 ) { cout << "Opening GLVis sockets." << endl; }
461 
462  socks_["A"] = new socketstream;
463  socks_["A"]->precision(8);
464 
465  socks_["B"] = new socketstream;
466  socks_["B"]->precision(8);
467 
468  socks_["H"] = new socketstream;
469  socks_["H"]->precision(8);
470 
471  if ( j_ )
472  {
473  socks_["J"] = new socketstream;
474  socks_["J"]->precision(8);
475  }
476  if ( k_ )
477  {
478  socks_["K"] = new socketstream;
479  socks_["K"]->precision(8);
480 
481  socks_["Psi"] = new socketstream;
482  socks_["Psi"]->precision(8);
483  }
484  if ( m_ )
485  {
486  socks_["M"] = new socketstream;
487  socks_["M"]->precision(8);
488  }
489  if ( myid_ == 0 ) { cout << "GLVis sockets open." << endl; }
490 }
491 
492 void
494 {
495  if (myid_ == 0) { cout << "Sending data to GLVis ..." << flush; }
496 
497  char vishost[] = "localhost";
498  int visport = 19916;
499 
500  int Wx = 0, Wy = 0; // window position
501  int Ww = 350, Wh = 350; // window size
502  int offx = Ww+10, offy = Wh+45; // window offsets
503 
504  VisualizeField(*socks_["A"], vishost, visport,
505  *a_, "Vector Potential (A)", Wx, Wy, Ww, Wh);
506  Wx += offx;
507 
508  VisualizeField(*socks_["B"], vishost, visport,
509  *b_, "Magnetic Flux Density (B)", Wx, Wy, Ww, Wh);
510  Wx += offx;
511 
512  VisualizeField(*socks_["H"], vishost, visport,
513  *h_, "Magnetic Field (H)", Wx, Wy, Ww, Wh);
514  Wx += offx;
515 
516  if ( j_ )
517  {
518  VisualizeField(*socks_["J"], vishost, visport,
519  *j_, "Current Density (J)", Wx, Wy, Ww, Wh);
520  }
521 
522  Wx = 0; Wy += offy; // next line
523 
524  if ( k_ )
525  {
526  VisualizeField(*socks_["K"], vishost, visport,
527  *k_, "Surface Current Density (K)", Wx, Wy, Ww, Wh);
528  Wx += offx;
529 
530  VisualizeField(*socks_["Psi"], vishost, visport,
531  *SurfCur_->GetPsi(),
532  "Surface Current Potential (Psi)", Wx, Wy, Ww, Wh);
533  Wx += offx;
534  }
535  if ( m_ )
536  {
537  VisualizeField(*socks_["M"], vishost, visport,
538  *m_, "Magnetization (M)", Wx, Wy, Ww, Wh);
539  // Wx += offx; // not used
540  }
541  if (myid_ == 0) { cout << " done." << endl; }
542 }
543 
546  Array<int> & kbcs,
547  Array<int> & vbcs, Vector & vbcv)
548  : H1FESpace_(&H1FESpace),
549  grad_(&grad),
550  kbcs_(&kbcs),
551  vbcs_(&vbcs),
552  vbcv_(&vbcv),
553  s0_(NULL),
554  psi_(NULL),
555  rhs_(NULL)
556 {
557  // Initialize MPI variables
558  MPI_Comm_rank(H1FESpace_->GetParMesh()->GetComm(), &myid_);
559 
560  s0_ = new ParBilinearForm(H1FESpace_);
561  s0_->AddBoundaryIntegrator(new DiffusionIntegrator);
562  s0_->Assemble();
563  s0_->Finalize();
564  S0_ = new HypreParMatrix;
565 
566  AttrToMarker(H1FESpace_->GetParMesh()->bdr_attributes.Max(),
567  *vbcs_, ess_bdr_);
568  H1FESpace_->GetEssentialTrueDofs(ess_bdr_, ess_bdr_tdofs_);
569 
570  non_k_bdr_.SetSize(H1FESpace_->GetParMesh()->bdr_attributes.Max());
571  non_k_bdr_ = 1;
572  for (int i=0; i<kbcs_->Size(); i++)
573  {
574  non_k_bdr_[(*kbcs_)[i]-1] = 0;
575  }
576 
577  psi_ = new ParGridFunction(H1FESpace_);
578  rhs_ = new ParGridFunction(H1FESpace_);
579 
580  pcg_ = NULL;
581  amg_ = NULL;
582 }
583 
585 {
586  delete psi_;
587  delete rhs_;
588 
589  delete pcg_;
590  delete amg_;
591 
592  delete S0_;
593 
594  delete s0_;
595 }
596 
597 void
599 {
600  delete pcg_;
601  delete amg_;
602 
603  amg_ = new HypreBoomerAMG(*S0_);
604  amg_->SetPrintLevel(0);
605  pcg_ = new HyprePCG(*S0_);
606  pcg_->SetTol(1e-14);
607  pcg_->SetMaxIter(200);
608  pcg_->SetPrintLevel(0);
609  pcg_->SetPreconditioner(*amg_);
610 }
611 
612 void
614 {
615  if (myid_ == 0) { cout << "Computing K ... " << flush; }
616 
617  // Apply piecewise constant voltage boundary condition
618  *psi_ = 0.0;
619  *rhs_ = 0.0;
620  Array<int> vbc_bdr_attr(H1FESpace_->GetParMesh()->bdr_attributes.Max());
621  for (int i=0; i<vbcs_->Size(); i++)
622  {
623  ConstantCoefficient voltage((*vbcv_)[i]);
624  vbc_bdr_attr = 0;
625  vbc_bdr_attr[(*vbcs_)[i]-1] = 1;
626  psi_->ProjectBdrCoefficient(voltage, vbc_bdr_attr);
627  }
628 
629  // Apply essential BC and form linear system
630  s0_->FormLinearSystem(ess_bdr_tdofs_, *psi_, *rhs_, *S0_, Psi_, RHS_);
631 
632  // Solve the linear system for Psi
633  if ( pcg_ == NULL ) { this->InitSolver(); }
634  pcg_->Mult(RHS_, Psi_);
635 
636  // Compute the parallel grid function corresponding to Psi
637  s0_->RecoverFEMSolution(Psi_, *rhs_, *psi_);
638 
639  // Compute the surface current from psi
640  grad_->Mult(*psi_, k);
641 
642  // Force the tangential part of k to be zero away from the intended surfaces
643  Vector vZero(3); vZero = 0.0;
644  VectorConstantCoefficient Zero(vZero);
645  k.ProjectBdrCoefficientTangent(Zero, non_k_bdr_);
646 
647  if (myid_ == 0) { cout << "done." << endl; }
648 }
649 
650 void
652 {
653  delete pcg_; pcg_ = NULL;
654  delete amg_; amg_ = NULL;
655  delete S0_; S0_ = new HypreParMatrix;
656 
657  psi_->Update();
658  rhs_->Update();
659 
660  s0_->Update();
661  s0_->Assemble();
662  s0_->Finalize();
663 
664  H1FESpace_->GetEssentialTrueDofs(ess_bdr_, ess_bdr_tdofs_);
665 }
666 
667 } // namespace electromagnetics
668 
669 } // namespace mfem
670 
671 #endif // MFEM_USE_MPI
const char vishost[]
void SetTol(double tol)
Definition: hypre.cpp:3990
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Definition: pfespace.cpp:1031
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:96
The Auxiliary-space Maxwell Solver in hypre.
Definition: hypre.hpp:1743
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
virtual int OrderW() const =0
Return the order of the determinant of the Jacobian (weight) of the transformation.
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:980
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
virtual void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
Form the linear system A X = B, corresponding to this bilinear form and the linear form b(...
void AttrToMarker(int max_attr, const Array< int > &attrs, Array< int > &marker)
Convert a set of attribute numbers to a marker array.
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve Ax=b with hypre&#39;s PCG.
Definition: hypre.cpp:4038
Integrator for (curl u, curl v) for Nedelec elements.
Vector coefficient that is constant in space and time.
int offx
virtual void Update(bool want_transform=true)
Definition: pfespace.cpp:3323
int Wx
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:68
virtual void SetIntRule(const IntegrationRule *ir)
Prescribe a fixed IntegrationRule to use (when ir != NULL) or let the integrator choose (when ir == N...
Definition: nonlininteg.hpp:57
Abstract parallel finite element space.
Definition: pfespace.hpp:28
void GetErrorEstimates(Vector &errors)
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: pgridfunc.cpp:561
STL namespace.
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
virtual void Mult(const Vector &x, Vector &y) const
Matrix multiplication: .
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:4010
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:480
virtual void Update(FiniteElementSpace *nfes=NULL)
Update the FiniteElementSpace and delete all data associated with the old one.
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds a domain integrator. Assumes ownership of bfi.
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1590
SurfaceCurrent(ParFiniteElementSpace &H1FESpace, ParDiscreteGradOperator &Grad, Array< int > &kbcs, Array< int > &vbcs, Vector &vbcv)
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition: pgridfunc.cpp:90
Jacobi preconditioner in hypre.
Definition: hypre.hpp:1421
void SetPrintLevel(int print_level)
Definition: hypre.hpp:1670
Data collection with VisIt I/O routines.
ParMesh * GetParMesh() const
Definition: pfespace.hpp:273
void Assemble(int skip_zeros=1)
void Assemble(int skip_zeros=1)
Assemble the local matrix.
const int visport
void SetTime(double t)
Set physical time (for time-dependent simulations)
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:281
void RegisterVisItFields(VisItDataCollection &visit_dc)
MPI_Comm GetComm() const
Definition: pmesh.hpp:351
void SetMaxIter(int max_iter)
Definition: hypre.cpp:4000
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition: fe_coll.hpp:380
A general vector function coefficient.
virtual void AddMult(const Vector &x, Vector &y, const double a=1.0) const
Add the matrix vector multiple to a vector: .
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:161
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:24
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:275
PCG solver in hypre.
Definition: hypre.hpp:1215
void ProjectBdrCoefficient(Coefficient *coeff[], VectorCoefficient *vcoeff, Array< int > &attr)
Definition: pgridfunc.cpp:669
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:41
virtual void Save() override
Save the collection and a VisIt root file.
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:687
HYPRE_Int HYPRE_BigInt
double L2ZZErrorEstimator(BilinearFormIntegrator &flux_integrator, const ParGridFunction &x, ParFiniteElementSpace &smooth_flux_fes, ParFiniteElementSpace &flux_fes, Vector &errors, int norm_p, double solver_tol, int solver_max_it)
Definition: pgridfunc.cpp:1270
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition: mesh.hpp:1023
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
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.
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:707
void SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
Definition: hypre.cpp:4015
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
void j_src(const Vector &x, double t, Vector &j)
Definition: maxwell.cpp:108
Class for parallel bilinear form.
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
void ComputeSurfaceCurrent(ParGridFunction &k)
virtual void AddMult(const Vector &x, Vector &y, const double a=1.0) const
Operator application: y+=A(x) (default) or y+=a*A(x).
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition: fe_coll.hpp:454
Vector data type.
Definition: vector.hpp:58
int offy
int Wy
int Ww
Class for parallel grid function.
Definition: pgridfunc.hpp:32
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
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
void SetSingularProblem()
Set this option when solving a curl-curl problem with zero mass term.
Definition: hypre.hpp:1802
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
virtual void RegisterField(const std::string &field_name, GridFunction *gf) override
Add a grid function to the collection and update the root file.
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
Definition: fespace.hpp:769