66 using namespace mfem::electromagnetics;
70 static double epsilon0_ = 8.8541878176e-12;
76 static Vector ds_params_(0);
81 static Vector cs_params_(0);
86 static Vector vp_params_(0);
91 static Vector e_uniform_(0);
97 int main(
int argc,
char *argv[])
104 const char *mesh_file =
"../../data/ball-nurbs.mesh";
107 int serial_ref_levels = 0;
108 int parallel_ref_levels = 0;
109 bool visualization =
true;
121 args.
AddOption(&mesh_file,
"-m",
"--mesh",
122 "Mesh file to use.");
124 "Finite element order (polynomial degree).");
125 args.
AddOption(&serial_ref_levels,
"-rs",
"--serial-ref-levels",
126 "Number of serial refinement levels.");
127 args.
AddOption(¶llel_ref_levels,
"-rp",
"--parallel-ref-levels",
128 "Number of parallel refinement levels.");
129 args.
AddOption(&e_uniform_,
"-uebc",
"--uniform-e-bc",
130 "Specify if the three components of the constant electric field");
131 args.
AddOption(&pw_eps_,
"-pwe",
"--piecewise-eps",
132 "Piecewise values of Permittivity");
133 args.
AddOption(&ds_params_,
"-ds",
"--dielectric-sphere-params",
134 "Center, Radius, and Permittivity of Dielectric Sphere");
135 args.
AddOption(&cs_params_,
"-cs",
"--charged-sphere-params",
136 "Center, Radius, and Total Charge of Charged Sphere");
137 args.
AddOption(&vp_params_,
"-vp",
"--voltaic-pile-params",
138 "Axis End Points, Radius, and Polarization of Cylindrical Voltaic Pile");
139 args.
AddOption(&dbcs,
"-dbcs",
"--dirichlet-bc-surf",
140 "Dirichlet Boundary Condition Surfaces");
141 args.
AddOption(&dbcv,
"-dbcv",
"--dirichlet-bc-vals",
142 "Dirichlet Boundary Condition Values");
143 args.
AddOption(&dbcg,
"-dbcg",
"--dirichlet-bc-gradient",
144 "-no-dbcg",
"--no-dirichlet-bc-gradient",
145 "Dirichlet Boundary Condition Gradient (phi = -z)");
146 args.
AddOption(&nbcs,
"-nbcs",
"--neumann-bc-surf",
147 "Neumann Boundary Condition Surfaces");
148 args.
AddOption(&nbcv,
"-nbcv",
"--neumann-bc-vals",
149 "Neumann Boundary Condition Values");
150 args.
AddOption(&maxit,
"-maxit",
"--max-amr-iterations",
151 "Max number of iterations in the main AMR loop.");
152 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
153 "--no-visualization",
154 "Enable or disable GLVis visualization.");
155 args.
AddOption(&visit,
"-visit",
"--visit",
"-no-visit",
"--no-visit",
156 "Enable or disable VisIt visualization.");
174 Mesh *mesh =
new Mesh(mesh_file, 1, 1);
179 cout <<
"Starting initialization." << endl;
186 if (serial_ref_levels > 0) { serial_ref_levels--; }
197 for (
int l = 0; l < serial_ref_levels; l++)
205 ParMesh pmesh(MPI_COMM_WORLD, *mesh);
209 int par_ref_levels = parallel_ref_levels;
210 for (
int l = 0; l < par_ref_levels; l++)
217 if ( dbcg && e_uniform_.
Size() != sdim )
221 e_uniform_(sdim-1) = 1.0;
225 if ( dbcv.
Size() < dbcs.
Size() && !dbcg )
242 VoltaSolver Volta(pmesh, order, dbcs, dbcv, nbcs, nbcv, *epsCoef,
260 if (mpi.
Root()) { cout <<
"Initialization done." << endl; }
267 const int max_dofs = 10000000;
268 for (
int it = 1; it <= maxit; it++)
272 cout <<
"\nAMR Iteration " << it << endl;
301 cout <<
"AMR iteration " << it <<
" complete." << endl;
305 if (prob_size > max_dofs)
309 cout <<
"Reached maximum number of dofs, exiting..." << endl;
316 if (mpi.
Root() && (it % 10 == 0))
318 cout <<
"press (q)uit or (c)ontinue --> " << flush;
321 MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
332 double local_max_err = errors.Max();
333 double global_max_err;
334 MPI_Allreduce(&local_max_err, &global_max_err, 1,
335 MPI_DOUBLE, MPI_MAX, pmesh.
GetComm());
339 const double frac = 0.7;
340 double threshold = frac * global_max_err;
341 if (mpi.
Root()) { cout <<
"Refining ..." << endl; }
349 if (mpi.
Root()) { cout <<
"Rebalancing ..." << endl; }
365 os <<
" ____ ____ __ __ " << endl
366 <<
" \\ \\ / /___ | |_/ |______ " << endl
367 <<
" \\ Y / _ \\| |\\ __\\__ \\ " << endl
368 <<
" \\ ( <_> ) |_| | / __ \\_ " << endl
369 <<
" \\___/ \\____/|____/__| (____ / " << endl
370 <<
" \\/ " << endl << flush;
380 if ( ds_params_.
Size() > 0 )
384 else if ( pw_eps_.
Size() > 0 )
403 for (
int i=0; i<x.
Size(); i++)
405 r2 += (x(i)-ds_params_(i))*(x(i)-ds_params_(i));
408 if ( sqrt(r2) <= ds_params_(x.
Size()) )
410 return ds_params_(x.
Size()+1) * epsilon0_;
423 if ( cs_params_(x.
Size()) > 0.0 )
428 rho = cs_params_(x.
Size()+1)/(M_PI*pow(cs_params_(x.
Size()),2));
431 rho = 0.75*cs_params_(x.
Size()+1)/(M_PI*pow(cs_params_(x.
Size()),3));
438 for (
int i=0; i<x.
Size(); i++)
440 r2 += (x(i)-cs_params_(i))*(x(i)-cs_params_(i));
443 if ( sqrt(r2) <= cs_params_(x.
Size()) )
463 for (
int i=0; i<x.
Size(); i++)
465 xu[i] -= vp_params_[i];
466 a[i] = vp_params_[x.
Size()+i] - vp_params_[i];
476 double r = vp_params_[2*x.
Size()];
486 if ( xa >= 0.0 && xa <= h*h && xp <= r )
488 p.
Add(vp_params_[2*x.
Size()+1]/h,a);
498 for (
int i=0; i<x.
Size(); i++)
500 phi -= x(i) * e_uniform_(i);
double phi_bc_uniform(const Vector &)
void WriteVisItFields(int it=0)
int Size() const
Logical size of the array.
Subclass constant coefficient.
void SetSize(int s)
Resize the vector if the new size is different.
HYPRE_Int GetProblemSize()
double Norml2() const
Returns the l2 norm of the vector.
int Size() const
Returns the size of the vector.
int GetNE() const
Returns number of elements.
void RegisterVisItFields(VisItDataCollection &visit_dc)
int main(int argc, char *argv[])
void voltaic_pile(const Vector &, Vector &)
void GetErrorEstimates(Vector &errors)
void Rebalance()
Load balance the mesh. NC meshes only.
bool Nonconforming() const
A simple convenience class that calls MPI_Init() at construction and MPI_Finalize() at destruction...
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Data collection with VisIt I/O routines.
void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
bool Root() const
Return true if WorldRank() == 0.
void PrintUsage(std::ostream &out) const
double charged_sphere(const Vector &)
int SpaceDimension() const
Base class Coefficient that may optionally depend on time.
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
bool RefineByError(const Array< double > &elem_error, double threshold, int nonconforming=-1, int nc_limit=0)
NURBSExtension * NURBSext
class for piecewise constant coefficient
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
void PrintOptions(std::ostream &out) const
void display_banner(ostream &os)
void EnsureNCMesh(bool triangles_nonconforming=false)
double dielectric_sphere(const Vector &)
class for C-function coefficient
Coefficient * SetupPermittivityCoefficient()
Class for parallel meshes.