62 using namespace mfem::electromagnetics;
65 static Vector ds_params_(0);
70 static Vector cs_params_(0);
75 static Vector vp_params_(0);
80 static Vector e_uniform_(0);
86 int main(
int argc,
char *argv[])
90 MPI_Init(&argc, &argv);
91 MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
92 MPI_Comm_rank(MPI_COMM_WORLD, &myid);
97 const char *mesh_file =
"butterfly_3d.mesh";
100 bool visualization =
true;
112 args.
AddOption(&mesh_file,
"-m",
"--mesh",
113 "Mesh file to use.");
115 "Finite element order (polynomial degree).");
116 args.
AddOption(&sr,
"-rs",
"--serial-ref-levels",
117 "Number of serial refinement levels.");
118 args.
AddOption(&pr,
"-rp",
"--parallel-ref-levels",
119 "Number of parallel refinement levels.");
120 args.
AddOption(&e_uniform_,
"-uebc",
"--uniform-e-bc",
121 "Specify if the three components of the constant electric field");
122 args.
AddOption(&ds_params_,
"-ds",
"--dielectric-sphere-params",
123 "Center, Radius, and Permittivity of Dielectric Sphere");
124 args.
AddOption(&cs_params_,
"-cs",
"--charged-sphere-params",
125 "Center, Radius, and Total Charge of Charged Sphere");
126 args.
AddOption(&vp_params_,
"-vp",
"--voltaic-pile-params",
127 "Axis End Points, Radius, and Polarization of Cylindrical Voltaic Pile");
128 args.
AddOption(&dbcs,
"-dbcs",
"--dirichlet-bc-surf",
129 "Dirichlet Boundary Condition Surfaces");
130 args.
AddOption(&dbcv,
"-dbcv",
"--dirichlet-bc-vals",
131 "Dirichlet Boundary Condition Values");
132 args.
AddOption(&dbcg,
"-dbcg",
"--dirichlet-bc-gradient",
133 "-no-dbcg",
"--no-dirichlet-bc-gradient",
134 "Dirichlet Boundary Condition Gradient (phi = -z)");
135 args.
AddOption(&nbcs,
"-nbcs",
"--neumann-bc-surf",
136 "Neumann Boundary Condition Surfaces");
137 args.
AddOption(&nbcv,
"-nbcv",
"--neumann-bc-vals",
138 "Neumann Boundary Condition Values");
139 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
140 "--no-visualization",
141 "Enable or disable GLVis visualization.");
142 args.
AddOption(&visit,
"-visit",
"--visit",
"-no-visit",
"--no-visit",
143 "Enable or disable VisIt visualization.");
163 ifstream imesh(mesh_file);
168 cerr <<
"\nCan not open mesh file: " << mesh_file <<
'\n' << endl;
173 mesh =
new Mesh(imesh, 1, 1);
181 if (myid == 0) { cout <<
"Starting initialization." << endl; }
184 if (mesh->
NURBSext && ref_levels < 2)
188 for (
int l = 0; l < ref_levels; l++)
205 ParMesh pmesh(MPI_COMM_WORLD, *mesh);
209 int par_ref_levels = pr;
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 )
239 VoltaSolver Volta(pmesh, order, dbcs, dbcv, nbcs, nbcv,
258 if (myid == 0) { cout <<
"Initialization done." << endl; }
265 const int max_dofs = 10000000;
266 for (
int it = 1; it <= 100; it++)
270 cout <<
"\nAMR Iteration " << it << endl;
293 if (myid == 0 && (visit || visualization)) { cout <<
"done." << endl; }
297 cout <<
"AMR iteration " << it <<
" complete." << endl;
301 if (prob_size > max_dofs)
305 cout <<
"Reached maximum number of dofs, exiting..." << endl;
312 if (myid == 0 && (it % 10 == 0))
314 cout <<
"press (q)uit or (c)ontinue --> " << flush;
317 MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
328 double local_max_err = errors.Max();
329 double global_max_err;
330 MPI_Allreduce(&local_max_err, &global_max_err, 1,
331 MPI_DOUBLE, MPI_MAX, pmesh.
GetComm());
336 const double frac = 0.7;
337 double threshold = frac * global_max_err;
338 for (
int i = 0; i < errors.Size(); i++)
340 if (errors[i] >= threshold) { ref_list.
Append(i); }
346 if (myid == 0) { cout <<
" Refinement ..." << flush; }
361 os <<
" ____ ____ __ __ " << endl
362 <<
" \\ \\ / /___ | |_/ |______ " << endl
363 <<
" \\ Y / _ \\| |\\ __\\__ \\ " << endl
364 <<
" \\ ( <_> ) |_| | / __ \\_ " << endl
365 <<
" \\___/ \\____/|____/__| (____ / " << endl
366 <<
" \\/ " << endl << flush;
376 for (
int i=0; i<x.
Size(); i++)
378 r2 += (x(i)-ds_params_(i))*(x(i)-ds_params_(i));
381 if ( sqrt(r2) <= ds_params_(x.
Size()) )
383 return ds_params_(x.
Size()+1) * epsilon0_;
396 if ( cs_params_(x.
Size()) > 0.0 )
401 rho = cs_params_(x.
Size()+1)/(M_PI*pow(cs_params_(x.
Size()),2));
404 rho = 0.75*cs_params_(x.
Size()+1)/(M_PI*pow(cs_params_(x.
Size()),3));
411 for (
int i=0; i<x.
Size(); i++)
413 r2 += (x(i)-cs_params_(i))*(x(i)-cs_params_(i));
416 if ( sqrt(r2) <= cs_params_(x.
Size()) )
436 for (
int i=0; i<x.
Size(); i++)
438 xu[i] -= vp_params_[i];
439 a[i] = vp_params_[x.
Size()+i] - vp_params_[i];
449 double r = vp_params_[2*x.
Size()];
459 if ( xa >= 0.0 && xa <= h*h && xp <= r )
461 p.
Add(vp_params_[2*x.
Size()+1]/h,a);
471 for (
int i=0; i<x.
Size(); i++)
473 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.
void SetSize(int s)
Resizes 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)
void voltaic_pile(const Vector &, Vector &)
void GetErrorEstimates(Vector &errors)
int Append(const T &el)
Append element to array, resize if necessary.
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)
void PrintUsage(std::ostream &out) const
double charged_sphere(const Vector &)
int SpaceDimension() const
int main(int argc, char *argv[])
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)
NURBSExtension * NURBSext
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
void PrintOptions(std::ostream &out) const
void display_banner(ostream &os)
double dielectric_sphere(const Vector &)
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Class for parallel meshes.