77 static Vector ds_params_(0);
82 static Vector cs_params_(0);
87 static Vector pc_params_(0);
90 static Vector vp_params_(0);
95 static Vector e_uniform_(0);
101 int main(
int argc,
char *argv[])
103 Mpi::Init(argc, argv);
109 const char *mesh_file =
"../../data/ball-nurbs.mesh";
112 int serial_ref_levels = 0;
113 int parallel_ref_levels = 0;
114 bool visualization =
true;
126 args.
AddOption(&mesh_file,
"-m",
"--mesh",
127 "Mesh file to use.");
129 "Finite element order (polynomial degree).");
130 args.
AddOption(&serial_ref_levels,
"-rs",
"--serial-ref-levels",
131 "Number of serial refinement levels.");
132 args.
AddOption(¶llel_ref_levels,
"-rp",
"--parallel-ref-levels",
133 "Number of parallel refinement levels.");
134 args.
AddOption(&e_uniform_,
"-uebc",
"--uniform-e-bc",
135 "Specify if the three components of the constant " 137 args.
AddOption(&pw_eps_,
"-pwe",
"--piecewise-eps",
138 "Piecewise values of Permittivity");
139 args.
AddOption(&ds_params_,
"-ds",
"--dielectric-sphere-params",
140 "Center, Radius, and Permittivity of Dielectric Sphere");
141 args.
AddOption(&cs_params_,
"-cs",
"--charged-sphere-params",
142 "Center, Radius, and Total Charge of Charged Sphere");
143 args.
AddOption(&pc_params_,
"-pc",
"--point-charge-params",
144 "Charges and locations of Point Charges");
145 args.
AddOption(&vp_params_,
"-vp",
"--voltaic-pile-params",
146 "Axis End Points, Radius, and " 147 "Polarization of Cylindrical Voltaic Pile");
148 args.
AddOption(&dbcs,
"-dbcs",
"--dirichlet-bc-surf",
149 "Dirichlet Boundary Condition Surfaces");
150 args.
AddOption(&dbcv,
"-dbcv",
"--dirichlet-bc-vals",
151 "Dirichlet Boundary Condition Values");
152 args.
AddOption(&dbcg,
"-dbcg",
"--dirichlet-bc-gradient",
153 "-no-dbcg",
"--no-dirichlet-bc-gradient",
154 "Dirichlet Boundary Condition Gradient (phi = -z)");
155 args.
AddOption(&nbcs,
"-nbcs",
"--neumann-bc-surf",
156 "Neumann Boundary Condition Surfaces");
157 args.
AddOption(&nbcv,
"-nbcv",
"--neumann-bc-vals",
158 "Neumann Boundary Condition Values");
159 args.
AddOption(&maxit,
"-maxit",
"--max-amr-iterations",
160 "Max number of iterations in the main AMR loop.");
161 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
162 "--no-visualization",
163 "Enable or disable GLVis visualization.");
164 args.
AddOption(&visit,
"-visit",
"--visit",
"-no-visit",
"--no-visit",
165 "Enable or disable VisIt visualization.");
183 Mesh *mesh =
new Mesh(mesh_file, 1, 1);
188 cout <<
"Starting initialization." << endl;
195 if (serial_ref_levels > 0) { serial_ref_levels--; }
206 for (
int l = 0; l < serial_ref_levels; l++)
214 ParMesh pmesh(MPI_COMM_WORLD, *mesh);
218 int par_ref_levels = parallel_ref_levels;
219 for (
int l = 0; l < par_ref_levels; l++)
228 if ( dbcg && e_uniform_.
Size() != sdim )
232 e_uniform_(sdim-1) = 1.0;
236 if ( dbcv.
Size() < dbcs.
Size() && !dbcg )
253 VoltaSolver Volta(pmesh, order, dbcs, dbcv, nbcs, nbcv, *epsCoef,
272 if (Mpi::Root()) { cout <<
"Initialization done." << endl; }
279 const int max_dofs = 10000000;
280 for (
int it = 1; it <= maxit; it++)
284 cout <<
"\nAMR Iteration " << it << endl;
313 cout <<
"AMR iteration " << it <<
" complete." << endl;
317 if (prob_size > max_dofs)
321 cout <<
"Reached maximum number of dofs, exiting..." << endl;
332 if (Mpi::Root() && (it % 10 == 0))
334 cout <<
"press (q)uit or (c)ontinue --> " << flush;
337 MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
348 double local_max_err = errors.Max();
349 double global_max_err;
350 MPI_Allreduce(&local_max_err, &global_max_err, 1,
351 MPI_DOUBLE, MPI_MAX, pmesh.
GetComm());
355 const double frac = 0.7;
356 double threshold = frac * global_max_err;
357 if (Mpi::Root()) { cout <<
"Refining ..." << endl; }
365 if (Mpi::Root()) { cout <<
"Rebalancing ..." << endl; }
381 os <<
" ____ ____ __ __ " << endl
382 <<
" \\ \\ / /___ | |_/ |______ " << endl
383 <<
" \\ Y / _ \\| |\\ __\\__ \\ " << endl
384 <<
" \\ ( <_> ) |_| | / __ \\_ " << endl
385 <<
" \\___/ \\____/|____/__| (____ / " << endl
386 <<
" \\/ " << endl << flush;
396 if ( ds_params_.
Size() > 0 )
400 else if ( pw_eps_.
Size() > 0 )
419 for (
int i=0; i<x.
Size(); i++)
421 r2 += (x(i)-ds_params_(i))*(x(i)-ds_params_(i));
424 if ( sqrt(r2) <= ds_params_(x.
Size()) )
426 return ds_params_(x.
Size()+1) * epsilon0_;
439 if ( cs_params_(x.
Size()) > 0.0 )
444 rho = cs_params_(x.
Size()+1) /
445 (M_PI * pow(cs_params_(x.
Size()), 2));
448 rho = 0.75 * cs_params_(x.
Size()+1) /
449 (M_PI * pow(cs_params_(x.
Size()), 3));
456 for (
int i=0; i<x.
Size(); i++)
458 r2 += (x(i) - cs_params_(i)) * (x(i) - cs_params_(i));
461 if ( sqrt(r2) <= cs_params_(x.
Size()) )
481 for (
int i=0; i<x.
Size(); i++)
483 xu[i] -= vp_params_[i];
484 a[i] = vp_params_[x.
Size()+i] - vp_params_[i];
487 double h =
a.Norml2();
494 double r = vp_params_[2 * x.
Size()];
499 xu.Add(-xa / (h * h),
a);
502 double xp = xu.Norml2();
504 if ( xa >= 0.0 && xa <= h*h && xp <= r )
506 p.Add(vp_params_[2 * x.
Size() + 1] / h,
a);
516 for (
int i=0; i<x.
Size(); i++)
518 phi -= x(i) * e_uniform_(i);
double phi_bc_uniform(const Vector &)
void WriteVisItFields(int it=0)
A coefficient that is constant across space and time.
void PrintOptions(std::ostream &out) const
Print the options.
void SetSize(int s)
Resize the vector to size s.
void PrintUsage(std::ostream &out) const
Print the usage message.
bool Nonconforming() const
int Size() const
Returns the size of the vector.
bool Good() const
Return true if the command line options were parsed successfully.
void RegisterVisItFields(VisItDataCollection &visit_dc)
void voltaic_pile(const Vector &, Vector &)
void GetErrorEstimates(Vector &errors)
void display_banner(ostream &os)
void Finalize(bool refine=false, bool fix_orientation=false) override
Finalize the construction of a general Mesh.
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void EnsureNCMesh(bool simplices_nonconforming=false)
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Data collection with VisIt I/O routines.
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Set the curvature of the mesh nodes using the given polynomial degree.
HYPRE_BigInt GetProblemSize()
double charged_sphere(const Vector &)
double p(const Vector &x, double t)
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
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)
Add a boolean option and set 'var' to receive the value. Enable/disable tags are used to set the bool...
int SpaceDimension() const
Dimension of the physical space containing the mesh.
bool RefineByError(const Array< double > &elem_error, double threshold, int nonconforming=-1, int nc_limit=0)
int GetNE() const
Returns number of elements.
NURBSExtension * NURBSext
Optional NURBS mesh extension.
A piecewise constant coefficient with the constants keyed off the element attribute numbers...
int Size() const
Return the logical size of the array.
double dielectric_sphere(const Vector &)
A general function coefficient.
int main(int argc, char *argv[])
Coefficient * SetupPermittivityCoefficient()
Class for parallel meshes.