22 namespace electromagnetics
25 VoltaSolver::VoltaSolver(
ParMesh & pmesh,
int order,
29 double (*phi_bc )(
const Vector&),
30 double (*rho_src)(
const Vector&),
70 rho_src_func_(rho_src),
72 point_charge_params_(point_charges),
76 MPI_Comm_size(pmesh_->
GetComm(), &num_procs_);
77 MPI_Comm_rank(pmesh_->
GetComm(), &myid_);
93 if ( phi_bc_func_ != NULL )
99 if ( rho_src_func_ != NULL )
105 if ( p_src_func_ != NULL )
139 if ( point_charge_params_.
Size() > 0 )
142 int npts = point_charge_params_.
Size() / (
dim + 1);
143 point_charges_.resize(npts);
146 for (
int i=0; i<npts; i++)
148 for (
int d=0; d<
dim; d++)
150 cent[d] = point_charge_params_[(
dim + 1) * i + d];
152 double s = point_charge_params_[(
dim + 1) * i +
dim];
155 point_charges_[i]->SetScale(
s);
156 point_charges_[i]->SetDeltaCenter(cent);
181 if ( nbcs_->
Size() > 0 )
214 delete hCurlHDivEps_;
219 delete HCurlFESpace_;
223 for (
unsigned int i=0; i<point_charges_.size(); i++)
225 delete point_charges_[i];
228 map<string,socketstream*>::iterator mit;
229 for (mit=socks_.begin(); mit!=socks_.end(); mit++)
250 cout <<
"Number of H1 unknowns: " << size_h1 << endl;
251 cout <<
"Number of H(Curl) unknowns: " << size_nd << endl;
252 cout <<
"Number of H(Div) unknowns: " << size_rt << endl;
253 cout <<
"Number of L2 unknowns: " << size_l2 << endl;
259 if (myid_ == 0) { cout <<
"Assembling ... " << flush; }
303 if (myid_ == 0) { cout <<
"done." << endl << flush; }
309 if (myid_ == 0) { cout <<
"Updating ..." << endl; }
314 H1FESpace_->
Update(
false);
315 HCurlFESpace_->
Update(
false);
316 HDivFESpace_->
Update(
false);
317 L2FESpace_->
Update(
false);
327 if ( rho_src_ ) { rho_src_->
Update(); }
328 if ( sigma_src_ ) { sigma_src_->
Update(); }
329 if ( p_src_ ) { p_src_->
Update(); }
336 if ( h1Mass_ ) { h1Mass_->
Update(); }
337 if ( h1SurfMass_ ) { h1SurfMass_->
Update(); }
338 if ( hCurlHDiv_ ) { hCurlHDiv_->
Update(); }
339 if ( weakDiv_ ) { weakDiv_->
Update(); }
349 if (myid_ == 0) { cout <<
"Running solver ... " << endl; }
354 if ( dbcs_->
Size() > 0 )
365 for (
int i=0; i<dbcs_->
Size(); i++)
369 if ((*dbcs_)[i] <= dbc_bdr_attr.Size())
371 dbc_bdr_attr[(*dbcs_)[i]-1] = 1;
382 h1Mass_->
AddMult(*rho_src_, *rhod_);
389 weakDiv_->
AddMult(*p_src_, *rhod_);
398 for (
int i=0; i<nbcs_->
Size(); i++)
402 if ((*nbcs_)[i] <= nbc_bdr_attr.Size())
404 nbc_bdr_attr[(*nbcs_)[i]-1] = 1;
408 h1SurfMass_->
AddMult(*sigma_src_, *rhod_);
412 if ( dbcs_->
Size() > 0 )
423 ess_bdr_tdofs_[0] = 0;
453 grad_->
Mult(*phi_, *e_); *e_ *= -1.0;
456 if (myid_ == 0) { cout <<
"Computing D ..." << flush; }
459 hCurlHDivEps_->
Mult(*e_, ed);
462 hCurlHDiv_->
AddMult(*p_src_, ed, -1.0);
482 div_->
Mult(*d_, *rho_);
484 if (myid_ == 0) { cout <<
"done." << flush; }
488 double charge_rho = (*l2_vol_int_)(*rho_);
491 double charge_D = (*rt_surf_int_)(*d_);
495 cout << endl <<
"Total charge: \n" 496 <<
" Volume integral of charge density: " << charge_rho
497 <<
"\n Surface integral of dielectric flux: " << charge_D
502 if (myid_ == 0) { cout <<
"Solver done. " << endl; }
508 if (myid_ == 0) { cout <<
"Estimating Error ... " << flush; }
522 smooth_flux_fes, flux_fes, errors, norm_p);
524 if (myid_ == 0) { cout <<
"done." << endl; }
530 visit_dc_ = &visit_dc;
536 if ( rho_src_ ) { visit_dc.
RegisterField(
"Rho Source", rho_src_); }
537 if ( p_src_ ) { visit_dc.
RegisterField(
"P Source", p_src_); }
538 if ( sigma_src_ ) { visit_dc.
RegisterField(
"Sigma Source", sigma_src_); }
546 if (myid_ == 0) { cout <<
"Writing VisIt files ..." << flush; }
553 if (myid_ == 0) { cout <<
" done." << endl; }
560 if ( myid_ == 0 ) { cout <<
"Opening GLVis sockets." << endl; }
563 socks_[
"Phi"]->precision(8);
566 socks_[
"D"]->precision(8);
569 socks_[
"E"]->precision(8);
572 socks_[
"Rho"]->precision(8);
577 socks_[
"RhoSrc"]->precision(8);
582 socks_[
"PSrc"]->precision(8);
587 socks_[
"SigmaSrc"]->precision(8);
594 if (myid_ == 0) { cout <<
"Sending data to GLVis ..." << flush; }
600 int Ww = 350,
Wh = 350;
604 *phi_,
"Electric Potential (Phi)",
Wx,
Wy,
Ww,
Wh);
608 *e_,
"Electric Field (E)",
Wx,
Wy,
Ww,
Wh);
612 *d_,
"Electric Displacement (D)",
Wx,
Wy,
Ww,
Wh);
616 *rho_,
"Charge Density",
Wx,
Wy,
Ww,
Wh);
622 *rho_src_,
"Charge Density Source (Rho)",
Wx,
Wy,
Ww,
Wh);
628 *p_src_,
"Electric Polarization Source (P)",
635 *sigma_src_,
"Surface Charge Density Source (Sigma)",
639 if (myid_ == 0) { cout <<
" done." << endl; }
646 #endif // MFEM_USE_MPI
void WriteVisItFields(int it=0)
Class for domain integration L(v) := (f, v)
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
A coefficient that is constant across space and time.
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.
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve Ax=b with hypre's PCG.
virtual void Update(bool want_transform=true)
int Size() const
Returns the size of the vector.
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Delta function coefficient optionally multiplied by a weight coefficient and a scaled time dependent ...
Abstract parallel finite element space.
void RegisterVisItFields(VisItDataCollection &visit_dc)
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
void SetPrintLevel(int print_lvl)
void GetErrorEstimates(Vector &errors)
The BoomerAMG solver in hypre.
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Jacobi preconditioner in hypre.
Data collection with VisIt I/O routines.
HYPRE_BigInt GetProblemSize()
void VisualizeField(socketstream &sock, const char *vishost, int visport, ParGridFunction &gf, const char *title, int x=0, int y=0, int w=400, int h=400, bool vec=false)
void SetTime(double t)
Set physical time (for time-dependent simulations)
HYPRE_BigInt GlobalTrueVSize() const
void SetMaxIter(int max_iter)
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
A general vector function coefficient.
Wrapper for hypre's parallel vector class.
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
void ProjectBdrCoefficient(Coefficient *coeff[], VectorCoefficient *vcoeff, Array< int > &attr)
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
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.
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)
int SpaceDimension() const
Dimension of the physical space containing the mesh.
void SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
int Size() const
Return the logical size of the array.
A general function coefficient.
Class for parallel grid function.
Wrapper for hypre's ParCSR matrix class.
virtual void Assemble(int skip_zeros=1)
Construct the internal matrix representation of the discrete linear operator.
Class for parallel meshes.
Arbitrary order "L2-conforming" discontinuous finite elements.
virtual void RegisterField(const std::string &field_name, GridFunction *gf) override
Add a grid function to the collection and update the root file.