21using namespace common;
23namespace electromagnetics
43 hDivHCurlMuInv_(NULL),
58 muInvCoef_(&muInvCoef),
67 MPI_Comm_size(pmesh_->
GetComm(), &num_procs_);
68 MPI_Comm_rank(pmesh_->
GetComm(), &myid_);
88 for (
int i=0; i<kbcs.
Size(); i++)
90 non_k_bdr_[kbcs[i]-1] = 0;
109 if ( j_src_ != NULL )
116 if ( m_src_ != NULL )
147 if ( jCoef_ || kbcs.
Size() > 0 )
156 irOrder, NULL, NULL, grad_);
159 if ( kbcs.
Size() > 0 )
200 delete curlMuInvCurl_;
202 delete hDivHCurlMuInv_;
203 delete weakCurlMuInv_;
206 delete HCurlFESpace_;
209 map<string,socketstream*>::iterator mit;
210 for (mit=socks_.begin(); mit!=socks_.end(); mit++)
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;
239 if (myid_ == 0) { cout <<
"Assembling ..." << flush; }
258 if ( weakCurlMuInv_ )
264 if (myid_ == 0) { cout <<
" done." << endl; }
270 if (myid_ == 0) { cout <<
"Updating ..." << endl; }
275 H1FESpace_->
Update(
false);
276 HCurlFESpace_->
Update(
false);
277 HDivFESpace_->
Update(
false);
287 if ( jr_ ) { jr_->
Update(); }
288 if ( j_ ) { j_->
Update(); }
289 if ( k_ ) { k_->
Update(); }
290 if ( m_ ) { m_->
Update(); }
295 hDivHCurlMuInv_->
Update();
296 if ( weakCurlMuInv_ ) { weakCurlMuInv_->
Update(); }
300 if ( grad_ ) { grad_->
Update(); }
301 if ( DivFreeProj_ ) { DivFreeProj_->
Update(); }
302 if ( SurfCur_ ) { SurfCur_->
Update(); }
308 if (myid_ == 0) { cout <<
"Running solver ... " << endl; }
332 DivFreeProj_->
Mult(*jr_, *j_);
335 hCurlMass_->
AddMult(*j_, *jd_);
342 weakCurlMuInv_->
AddMult(*m_, *jd_, mu0_);
356 HypreAMS ams(CurlMuInvCurl, HCurlFESpace_);
374 curl_->
Mult(*a_, *b_);
377 if (myid_ == 0) { cout <<
"Computing H ... " << flush; }
379 hDivHCurlMuInv_->
Mult(*b_, *bd_);
382 hDivHCurlMuInv_->
AddMult(*m_, *bd_, -1.0 * mu0_);
401 if (myid_ == 0) { cout <<
"done." << flush; }
403 if (myid_ == 0) { cout <<
" Solver done. " << endl; }
409 if (myid_ == 0) { cout <<
"Estimating Error ... " << flush; }
422 smooth_flux_fes, flux_fes, errors, norm_p);
424 if (myid_ == 0) { cout <<
"done." << endl; }
430 visit_dc_ = &visit_dc;
446 if (myid_ == 0) { cout <<
"Writing VisIt files ..." << flush; }
453 if (myid_ == 0) { cout <<
" done." << endl; }
460 if ( myid_ == 0 ) { cout <<
"Opening GLVis sockets." << endl; }
463 socks_[
"A"]->precision(8);
466 socks_[
"B"]->precision(8);
469 socks_[
"H"]->precision(8);
474 socks_[
"J"]->precision(8);
479 socks_[
"K"]->precision(8);
482 socks_[
"Psi"]->precision(8);
487 socks_[
"M"]->precision(8);
489 if ( myid_ == 0 ) { cout <<
"GLVis sockets open." << endl; }
495 if (myid_ == 0) { cout <<
"Sending data to GLVis ..." << flush; }
500 int Ww = 350,
Wh = 350;
504 *a_,
"Vector Potential (A)",
Wx,
Wy,
Ww,
Wh);
508 *b_,
"Magnetic Flux Density (B)",
Wx,
Wy,
Ww,
Wh);
512 *h_,
"Magnetic Field (H)",
Wx,
Wy,
Ww,
Wh);
518 *j_,
"Current Density (J)",
Wx,
Wy,
Ww,
Wh);
526 *k_,
"Surface Current Density (K)",
Wx,
Wy,
Ww,
Wh);
531 "Surface Current Potential (Psi)",
Wx,
Wy,
Ww,
Wh);
537 *m_,
"Magnetization (M)",
Wx,
Wy,
Ww,
Wh);
540 if (myid_ == 0) { cout <<
" done." << endl; }
547 : H1FESpace_(&H1FESpace),
571 for (
int i=0; i<kbcs_->
Size(); i++)
573 non_k_bdr_[(*kbcs_)[i]-1] = 0;
614 if (myid_ == 0) { cout <<
"Computing K ... " << flush; }
620 for (
int i=0; i<vbcs_->
Size(); i++)
624 vbc_bdr_attr[(*vbcs_)[i]-1] = 1;
633 pcg_->
Mult(RHS_, Psi_);
639 grad_->
Mult(*psi_, k);
642 Vector vZero(3); vZero = 0.0;
646 if (myid_ == 0) { cout <<
"done." << endl; }
652 delete pcg_; pcg_ = NULL;
653 delete amg_; amg_ = NULL;
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
int Size() const
Return the logical size of the array.
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
A coefficient that is constant across space and time.
Integrator for for Nedelec elements.
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
void SetTime(real_t t)
Set physical time (for time-dependent simulations)
virtual void Assemble(int skip_zeros=1)
Construct the internal matrix representation of the discrete linear operator.
The Auxiliary-space Maxwell Solver in hypre.
void SetSingularProblem()
Set this option when solving a curl-curl problem with zero mass term.
The BoomerAMG solver in hypre.
void SetPrintLevel(int print_level)
Jacobi preconditioner in hypre.
void Mult(const HypreParVector &b, HypreParVector &x) const override
Solve Ax=b with hypre's PCG.
void SetPrintLevel(int print_lvl)
void SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
void SetMaxIter(int max_iter)
Wrapper for hypre's ParCSR matrix class.
Wrapper for hypre's parallel vector class.
Class for an integration rule - an Array of IntegrationPoint.
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
virtual void SetIntRule(const IntegrationRule *ir)
Prescribe a fixed IntegrationRule to use, or set to null to let the integrator choose an appropriate ...
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Geometry::Type GetTypicalElementGeometry() const
If the local mesh is not empty, return GetElementGeometry(0); otherwise, return a typical Geometry pr...
int Dimension() const
Dimension of the reference space used within the elements.
ElementTransformation * GetTypicalElementTransformation()
If the local mesh is not empty return GetElementTransformation(0); otherwise, return the identity tra...
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Arbitrary order H(curl)-conforming Nedelec finite elements.
Abstract parallel finite element space.
void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1) const override
HYPRE_BigInt GlobalTrueVSize() const
ParMesh * GetParMesh() const
void Update(bool want_transform=true) override
Class for parallel grid function.
void ProjectBdrCoefficient(Coefficient *coeff[], VectorCoefficient *vcoeff, const Array< int > &attr)
void ProjectCoefficient(Coefficient &coeff) override
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
void Update() override
Transform by the Space UpdateMatrix (e.g., on Mesh change).
void ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff, const Array< int > &bdr_attr) override
Project the tangential components of the given VectorCoefficient on the boundary. Only boundary attri...
Class for parallel meshes.
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Vector coefficient that is constant in space and time.
A general vector function coefficient.
Data collection with VisIt I/O routines.
void Save() override
Save the collection and a VisIt root file.
void RegisterField(const std::string &field_name, GridFunction *gf) override
Add a grid function to the collection and update the root file.
void Mult(const Vector &x, Vector &y) const override
Operator application: y=A(x).
void ComputeSurfaceCurrent(ParGridFunction &k)
SurfaceCurrent(ParFiniteElementSpace &H1FESpace, ParDiscreteGradOperator &Grad, Array< int > &kbcs, Array< int > &vbcs, Vector &vbcv)
ParGridFunction * GetPsi()
void GetErrorEstimates(Vector &errors)
HYPRE_BigInt GetProblemSize()
void WriteVisItFields(int it=0)
void DisplayToGLVis(int visport=19916)
TeslaSolver(ParMesh &pmesh, int order, Array< int > &kbcs, Array< int > &vbcs, Vector &vbcv, Coefficient &muInvCoef, void(*a_bc)(const Vector &, Vector &), void(*j_src)(const Vector &, Vector &), void(*m_src)(const Vector &, Vector &))
void RegisterVisItFields(VisItDataCollection &visit_dc)
void j_src(const Vector &x, real_t t, Vector &j)
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)
void AttrToMarker(int max_attr, const Array< int > &attrs, Array< int > &marker)
Convert a set of attribute numbers to a marker array.
real_t L2ZZErrorEstimator(BilinearFormIntegrator &flux_integrator, const ParGridFunction &x, ParFiniteElementSpace &smooth_flux_fes, ParFiniteElementSpace &flux_fes, Vector &errors, int norm_p, real_t solver_tol, int solver_max_it)
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)