68 using namespace mfem::electromagnetics;
74 static Vector pw_mu_inv_(0);
75 static Vector ms_params_(0);
81 static Vector cr_params_(0);
87 static Vector bm_params_(0);
91 static Vector ha_params_(0);
98 static Vector b_uniform_(0);
107 int main(
int argc,
char *argv[])
114 const char *mesh_file =
"../../data/ball-nurbs.mesh";
117 int serial_ref_levels = 0;
118 int parallel_ref_levels = 0;
119 bool visualization =
true;
128 args.
AddOption(&mesh_file,
"-m",
"--mesh",
129 "Mesh file to use.");
131 "Finite element order (polynomial degree).");
132 args.
AddOption(&serial_ref_levels,
"-rs",
"--serial-ref-levels",
133 "Number of serial refinement levels.");
134 args.
AddOption(¶llel_ref_levels,
"-rp",
"--parallel-ref-levels",
135 "Number of parallel refinement levels.");
136 args.
AddOption(&b_uniform_,
"-ubbc",
"--uniform-b-bc",
137 "Specify if the three components of the constant magnetic flux density");
138 args.
AddOption(&pw_mu_,
"-pwm",
"--piecewise-mu",
139 "Piecewise values of Permeability");
140 args.
AddOption(&ms_params_,
"-ms",
"--magnetic-shell-params",
141 "Center, Inner Radius, Outer Radius, and Permeability of Magnetic Shell");
142 args.
AddOption(&cr_params_,
"-cr",
"--current-ring-params",
143 "Axis End Points, Inner Radius, Outer Radius and Total Current of Annulus");
144 args.
AddOption(&bm_params_,
"-bm",
"--bar-magnet-params",
145 "Axis End Points, Radius, and Magnetic Field of Cylindrical Magnet");
146 args.
AddOption(&ha_params_,
"-ha",
"--halbach-array-params",
147 "Bounding Box Corners and Number of Segments");
148 args.
AddOption(&kbcs,
"-kbcs",
"--surface-current-bc",
149 "Surfaces for the Surface Current (K) Boundary Condition");
150 args.
AddOption(&vbcs,
"-vbcs",
"--voltage-bc-surf",
151 "Voltage Boundary Condition Surfaces (to drive K)");
152 args.
AddOption(&vbcv,
"-vbcv",
"--voltage-bc-vals",
153 "Voltage Boundary Condition Values (to drive K)");
154 args.
AddOption(&maxit,
"-maxit",
"--max-amr-iterations",
155 "Max number of iterations in the main AMR loop.");
156 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
157 "--no-visualization",
158 "Enable or disable GLVis visualization.");
159 args.
AddOption(&visit,
"-visit",
"--visit",
"-no-visit",
"--no-visit",
160 "Enable or disable VisIt visualization.");
178 Mesh *mesh =
new Mesh(mesh_file, 1, 1);
182 cout <<
"Starting initialization." << endl;
189 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++)
220 if ( ( vbcs.
Size() > 0 && kbcs.
Size() == 0 ) ||
221 ( kbcs.
Size() > 0 && vbcs.
Size() == 0 ) ||
226 cout <<
"The surface current (K) boundary condition requires "
227 <<
"surface current boundary condition surfaces (with -kbcs), "
228 <<
"voltage boundary condition surface (with -vbcs), "
229 <<
"and voltage boundary condition values (with -vbcv)."
239 TeslaSolver Tesla(pmesh, order, kbcs, vbcs, vbcv, *muInvCoef,
258 if (mpi.
Root()) { cout <<
"Initialization done." << endl; }
265 const int max_dofs = 10000000;
266 for (
int it = 1; it <= maxit; it++)
270 cout <<
"\nAMR Iteration " << it << endl;
299 cout <<
"AMR iteration " << it <<
" complete." << endl;
303 if (prob_size > max_dofs)
307 cout <<
"Reached maximum number of dofs, exiting..." << endl;
318 if (mpi.
Root() && (it % 10 == 0))
320 cout <<
"press (q)uit or (c)ontinue --> " << flush;
323 MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
334 double local_max_err = errors.Max();
335 double global_max_err;
336 MPI_Allreduce(&local_max_err, &global_max_err, 1,
337 MPI_DOUBLE, MPI_MAX, pmesh.
GetComm());
341 const double frac = 0.5;
342 double threshold = frac * global_max_err;
343 if (mpi.
Root()) { cout <<
"Refining ..." << endl; }
351 if (mpi.
Root()) { cout <<
"Rebalancing ..." << endl; }
367 os <<
" ___________ __ " << endl
368 <<
" \\__ ___/___ _____| | _____ " << endl
369 <<
" | |_/ __ \\ / ___/ | \\__ \\ " << endl
370 <<
" | |\\ ___/ \\___ \\| |__/ __ \\_ " << endl
371 <<
" |____| \\___ >____ >____(____ / " << endl
372 <<
" \\/ \\/ \\/ " << endl << flush;
382 if ( ms_params_.
Size() > 0 )
386 else if ( pw_mu_.
Size() > 0 )
389 for (
int i = 0; i < pw_mu_.
Size(); i++)
391 MFEM_ASSERT( pw_mu_[i] > 0.0,
"permeability values must be positive" );
392 pw_mu_inv_[i] = 1.0/pw_mu_[i];
411 for (
int i = 0; i < x.
Size(); i++)
413 r2 += (x(i) - ms_params_(i))*(x(i) - ms_params_(i));
416 if ( sqrt(r2) >= ms_params_(x.
Size()) &&
417 sqrt(r2) <= ms_params_(x.
Size()+1) )
419 return mu0_*ms_params_(x.
Size()+2);
428 MFEM_ASSERT(x.
Size() == 3,
"current_ring source requires 3D space.");
439 for (
int i=0; i<x.
Size(); i++)
441 xu[i] -= cr_params_[i];
442 a[i] = cr_params_[x.
Size()+i] - cr_params_[i];
452 double ra = cr_params_[2*x.
Size()+0];
453 double rb = cr_params_[2*x.
Size()+1];
469 if ( xa >= 0.0 && xa <= h*h && xp >= ra && xp <= rb )
471 ju(0) = a(1) * xu(2) - a(2) * xu(1);
472 ju(1) = a(2) * xu(0) - a(0) * xu(2);
473 ju(2) = a(0) * xu(1) - a(1) * xu(0);
476 j.
Add(cr_params_[2*x.
Size()+2]/(h*(rb-ra)),ju);
493 for (
int i=0; i<x.
Size(); i++)
495 xu[i] -= bm_params_[i];
496 a[i] = bm_params_[x.
Size()+i] - bm_params_[i];
506 double r = bm_params_[2*x.
Size()];
516 if ( xa >= 0.0 && xa <= h*h && xp <= r )
518 m.
Add(bm_params_[2*x.
Size()+1]/h,a);
531 if ( x[0] < ha_params_[0] || x[0] > ha_params_[3] ||
532 x[1] < ha_params_[1] || x[1] > ha_params_[4] ||
533 x[2] < ha_params_[2] || x[2] > ha_params_[5] )
538 int ai = (int)ha_params_[6];
539 int ri = (int)ha_params_[7];
540 int n = (int)ha_params_[8];
542 int i = (int)n * (x[ai] - ha_params_[ai]) /
543 (ha_params_[ai+3] - ha_params_[ai]);
545 m[(ri + 1 + (i % 2)) % 3] = pow(-1.0,i/2);
553 a(0) = b_uniform_(1) * x(2);
554 a(1) = b_uniform_(2) * x(0);
555 a(2) = b_uniform_(0) * x(1);
562 return -x(x.
Size()-1);
int WorldSize() const
Return MPI_COMM_WORLD's size.
int Size() const
Logical size of the array.
Subclass constant coefficient.
void SetSize(int s)
Resize the vector to size s.
double Norml2() const
Returns the l2 norm of the vector.
int Size() const
Returns the size of the vector.
void bar_magnet(const Vector &, Vector &)
int GetNE() const
Returns number of elements.
void GetErrorEstimates(Vector &errors)
int main(int argc, char *argv[])
double phi_m_bc_uniform(const Vector &x)
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...
double magnetic_shell_inv(const Vector &x)
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
void current_ring(const Vector &, Vector &)
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.
HYPRE_Int GetProblemSize()
void PrintUsage(std::ostream &out) const
void RegisterVisItFields(VisItDataCollection &visit_dc)
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)
void a_bc_uniform(const Vector &, Vector &)
NURBSExtension * NURBSext
Optional NURBS mesh extension.
class for piecewise constant coefficient
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
void WriteVisItFields(int it=0)
Coefficient * SetupInvPermeabilityCoefficient()
void PrintOptions(std::ostream &out) const
double magnetic_shell(const Vector &)
virtual void Finalize(bool refine=false, bool fix_orientation=false)
Finalize the construction of a general Mesh.
void EnsureNCMesh(bool triangles_nonconforming=false)
class for C-function coefficient
void display_banner(ostream &os)
Class for parallel meshes.
void halbach_array(const Vector &, Vector &)