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 "
132 args.
AddOption(&pw_eps_,
"-pwe",
"--piecewise-eps",
133 "Piecewise values of Permittivity");
134 args.
AddOption(&ds_params_,
"-ds",
"--dielectric-sphere-params",
135 "Center, Radius, and Permittivity of Dielectric Sphere");
136 args.
AddOption(&cs_params_,
"-cs",
"--charged-sphere-params",
137 "Center, Radius, and Total Charge of Charged Sphere");
138 args.
AddOption(&vp_params_,
"-vp",
"--voltaic-pile-params",
139 "Axis End Points, Radius, and "
140 "Polarization of Cylindrical Voltaic Pile");
141 args.
AddOption(&dbcs,
"-dbcs",
"--dirichlet-bc-surf",
142 "Dirichlet Boundary Condition Surfaces");
143 args.
AddOption(&dbcv,
"-dbcv",
"--dirichlet-bc-vals",
144 "Dirichlet Boundary Condition Values");
145 args.
AddOption(&dbcg,
"-dbcg",
"--dirichlet-bc-gradient",
146 "-no-dbcg",
"--no-dirichlet-bc-gradient",
147 "Dirichlet Boundary Condition Gradient (phi = -z)");
148 args.
AddOption(&nbcs,
"-nbcs",
"--neumann-bc-surf",
149 "Neumann Boundary Condition Surfaces");
150 args.
AddOption(&nbcv,
"-nbcv",
"--neumann-bc-vals",
151 "Neumann Boundary Condition Values");
152 args.
AddOption(&maxit,
"-maxit",
"--max-amr-iterations",
153 "Max number of iterations in the main AMR loop.");
154 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
155 "--no-visualization",
156 "Enable or disable GLVis visualization.");
157 args.
AddOption(&visit,
"-visit",
"--visit",
"-no-visit",
"--no-visit",
158 "Enable or disable VisIt visualization.");
176 Mesh *mesh =
new Mesh(mesh_file, 1, 1);
181 cout <<
"Starting initialization." << endl;
188 if (serial_ref_levels > 0) { serial_ref_levels--; }
199 for (
int l = 0; l < serial_ref_levels; l++)
207 ParMesh pmesh(MPI_COMM_WORLD, *mesh);
211 int par_ref_levels = parallel_ref_levels;
212 for (
int l = 0; l < par_ref_levels; l++)
219 if ( dbcg && e_uniform_.
Size() != sdim )
223 e_uniform_(sdim-1) = 1.0;
227 if ( dbcv.
Size() < dbcs.
Size() && !dbcg )
244 VoltaSolver Volta(pmesh, order, dbcs, dbcv, nbcs, nbcv, *epsCoef,
262 if (mpi.
Root()) { cout <<
"Initialization done." << endl; }
269 const int max_dofs = 10000000;
270 for (
int it = 1; it <= maxit; it++)
274 cout <<
"\nAMR Iteration " << it << endl;
303 cout <<
"AMR iteration " << it <<
" complete." << endl;
307 if (prob_size > max_dofs)
311 cout <<
"Reached maximum number of dofs, exiting..." << endl;
322 if (mpi.
Root() && (it % 10 == 0))
324 cout <<
"press (q)uit or (c)ontinue --> " << flush;
327 MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
338 double local_max_err = errors.Max();
339 double global_max_err;
340 MPI_Allreduce(&local_max_err, &global_max_err, 1,
341 MPI_DOUBLE, MPI_MAX, pmesh.
GetComm());
345 const double frac = 0.7;
346 double threshold = frac * global_max_err;
347 if (mpi.
Root()) { cout <<
"Refining ..." << endl; }
355 if (mpi.
Root()) { cout <<
"Rebalancing ..." << endl; }
371 os <<
" ____ ____ __ __ " << endl
372 <<
" \\ \\ / /___ | |_/ |______ " << endl
373 <<
" \\ Y / _ \\| |\\ __\\__ \\ " << endl
374 <<
" \\ ( <_> ) |_| | / __ \\_ " << endl
375 <<
" \\___/ \\____/|____/__| (____ / " << endl
376 <<
" \\/ " << endl << flush;
386 if ( ds_params_.
Size() > 0 )
390 else if ( pw_eps_.
Size() > 0 )
409 for (
int i=0; i<x.
Size(); i++)
411 r2 += (x(i)-ds_params_(i))*(x(i)-ds_params_(i));
414 if ( sqrt(r2) <= ds_params_(x.
Size()) )
416 return ds_params_(x.
Size()+1) * epsilon0_;
429 if ( cs_params_(x.
Size()) > 0.0 )
434 rho = cs_params_(x.
Size()+1) /
435 (M_PI * pow(cs_params_(x.
Size()), 2));
438 rho = 0.75 * cs_params_(x.
Size()+1) /
439 (M_PI * pow(cs_params_(x.
Size()), 3));
446 for (
int i=0; i<x.
Size(); i++)
448 r2 += (x(i) - cs_params_(i)) * (x(i) - cs_params_(i));
451 if ( sqrt(r2) <= cs_params_(x.
Size()) )
471 for (
int i=0; i<x.
Size(); i++)
473 xu[i] -= vp_params_[i];
474 a[i] = vp_params_[x.
Size()+i] - vp_params_[i];
484 double r = vp_params_[2 * x.
Size()];
489 xu.
Add(-xa / (h * h), a);
494 if ( xa >= 0.0 && xa <= h*h && xp <= r )
496 p.
Add(vp_params_[2 * x.
Size() + 1] / h, a);
506 for (
int i=0; i<x.
Size(); i++)
508 phi -= x(i) * e_uniform_(i);
double phi_bc_uniform(const Vector &)
int WorldSize() const
Return MPI_COMM_WORLD's size.
void WriteVisItFields(int it=0)
int Size() const
Logical size of the array.
Subclass constant coefficient.
void SetSize(int s)
Resize the vector to size s.
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
Optional NURBS mesh extension.
class for piecewise constant coefficient
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
void PrintOptions(std::ostream &out) const
void EnsureNCMesh(bool triangles_nonconforming=false)
double dielectric_sphere(const Vector &)
class for C-function coefficient
Coefficient * SetupPermittivityCoefficient()
void display_banner(ostream &os)
Class for parallel meshes.