45 int main(
int argc,
char *argv[])
48 Mpi::Init(argc, argv);
49 int num_procs = Mpi::WorldSize();
50 int myid = Mpi::WorldRank();
54 const char *mesh_file =
"../data/inline-quad.mesh";
55 int ser_ref_levels = 2;
56 int par_ref_levels = 1;
59 bool visualization = 1;
62 args.
AddOption(&mesh_file,
"-m",
"--mesh",
64 args.
AddOption(&ser_ref_levels,
"-rs",
"--refine-serial",
65 "Number of times to refine the mesh uniformly in serial.");
66 args.
AddOption(&par_ref_levels,
"-rp",
"--refine-parallel",
67 "Number of times to refine the mesh uniformly in parallel.");
69 "Finite element order (polynomial degree).");
70 args.
AddOption(&
freq,
"-f",
"--frequency",
"Set the frequency for the exact"
72 args.
AddOption(&use_ams,
"-ams",
"--hypre-ams",
"-slu",
73 "--superlu",
"Use AMS or SuperLU solver.");
74 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
76 "Enable or disable GLVis visualization.");
84 Mesh *mesh =
new Mesh(mesh_file, 1, 1);
90 for (
int lev = 0; lev < ser_ref_levels; lev++)
99 ParMesh pmesh(MPI_COMM_WORLD, *mesh);
101 for (
int lev = 0; lev < par_ref_levels; lev++)
123 if (Mpi::Root()) { cout <<
"Number of H(Curl) unknowns: " << size << endl; }
160 sigmaMat(0,0) = 2.0; sigmaMat(1,1) = 2.0; sigmaMat(2,2) = 2.0;
161 sigmaMat(0,2) = 0.0; sigmaMat(2,0) = 0.0;
162 sigmaMat(0,1) = M_SQRT1_2; sigmaMat(1,0) = M_SQRT1_2;
163 sigmaMat(1,2) = M_SQRT1_2; sigmaMat(2,1) = M_SQRT1_2;
187 cout <<
"Size of linear system: "
201 #ifdef MFEM_USE_SUPERLU
205 cout <<
"Size of linear system: "
216 if (Mpi::Root()) { cout <<
"No solvers available." << endl; }
230 cout <<
"\n|| E_h - E ||_{H(Curl)} = " << error <<
'\n' << endl;
238 ostringstream mesh_name, sol_name;
239 mesh_name <<
"mesh." << setfill(
'0') << setw(6) << myid;
240 sol_name <<
"sol." << setfill(
'0') << setw(6) << myid;
242 ofstream mesh_ofs(mesh_name.str().c_str());
243 mesh_ofs.precision(8);
244 pmesh.
Print(mesh_ofs);
246 ofstream sol_ofs(sol_name.str().c_str());
247 sol_ofs.precision(8);
270 dy_sock.precision(8);
271 dz_sock.precision(8);
273 Vector xVec(3); xVec = 0.0; xVec(0) = 1;
274 Vector yVec(3); yVec = 0.0; yVec(1) = 1;
275 Vector zVec(3); zVec = 0.0; zVec(2) = 1;
301 x_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n"
302 <<
"solution\n" << pmesh << xComp << flush
303 <<
"window_title 'X component'" << endl;
304 y_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n"
305 <<
"solution\n" << pmesh << yComp << flush
306 <<
"window_geometry 403 0 400 350 "
307 <<
"window_title 'Y component'" << endl;
308 z_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n"
309 <<
"solution\n" << pmesh << zComp << flush
310 <<
"window_geometry 806 0 400 350 "
311 <<
"window_title 'Z component'" << endl;
319 dy_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n"
320 <<
"solution\n" << pmesh << dyComp << flush
321 <<
"window_geometry 403 375 400 350 "
322 <<
"window_title 'Y component of Curl'" << endl;
323 dz_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n"
324 <<
"solution\n" << pmesh << dzComp << flush
325 <<
"window_geometry 806 375 400 350 "
326 <<
"window_title 'Z component of Curl'" << endl;
336 xyMat(0,0) = 1.0; xyMat(1,1) = 1.0;
338 Vector zVec(3); zVec = 0.0; zVec(2) = 1;
363 xy_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n";
364 xy_sock.precision(8);
365 xy_sock <<
"solution\n" << pmesh << xyComp
366 <<
"window_title 'XY components'\n" << flush;
367 z_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n"
368 <<
"solution\n" << pmesh << zComp << flush
369 <<
"window_geometry 403 0 400 350 "
370 <<
"window_title 'Z component'" << endl;
378 dxy_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n"
379 <<
"solution\n" << pmesh << dxyComp << flush
380 <<
"window_geometry 0 375 400 350 "
381 <<
"window_title 'XY components of Curl'" << endl;
382 dz_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n"
383 <<
"solution\n" << pmesh << dzComp << flush
384 <<
"window_geometry 403 375 400 350 "
385 <<
"window_title 'Z component of Curl'" << endl;
400 sol_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n";
401 sol_sock.precision(8);
402 sol_sock <<
"solution\n" << pmesh << sol
403 <<
"window_title 'Solution'" << flush << endl;
404 dsol_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n"
405 <<
"solution\n" << pmesh << dsol << flush
406 <<
"window_geometry 0 375 400 350 "
407 <<
"window_title 'Curl of solution'" << endl;
422 E(0) = 1.1 * sin(
kappa * x(0) + 0.0 * M_PI);
423 E(1) = 1.2 * sin(
kappa * x(0) + 0.4 * M_PI);
424 E(2) = 1.3 * sin(
kappa * x(0) + 0.9 * M_PI);
428 E(0) = 1.1 * sin(
kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.0 * M_PI);
429 E(1) = 1.2 * sin(
kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.4 * M_PI);
430 E(2) = 1.3 * sin(
kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.9 * M_PI);
434 E(0) = 1.1 * sin(
kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.0 * M_PI);
435 E(1) = 1.2 * sin(
kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.4 * M_PI);
436 E(2) = 1.3 * sin(
kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.9 * M_PI);
437 E *= cos(
kappa * x(2));
445 double c4 = cos(
kappa * x(0) + 0.4 * M_PI);
446 double c9 = cos(
kappa * x(0) + 0.9 * M_PI);
455 double c0 = cos(
kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.0 * M_PI);
456 double c4 = cos(
kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.4 * M_PI);
457 double c9 = cos(
kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.9 * M_PI);
461 dE(2) = 1.2 * c4 - 1.1 * c0;
462 dE *=
kappa * M_SQRT1_2;
466 double s0 = sin(
kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.0 * M_PI);
467 double c0 = cos(
kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.0 * M_PI);
468 double s4 = sin(
kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.4 * M_PI);
469 double c4 = cos(
kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.4 * M_PI);
470 double c9 = cos(
kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.9 * M_PI);
471 double sk = sin(
kappa * x(2));
472 double ck = cos(
kappa * x(2));
474 dE(0) = 1.2 * s4 * sk + 1.3 * M_SQRT1_2 * c9 * ck;
475 dE(1) = -1.1 * s0 * sk - 1.3 * M_SQRT1_2 * c9 * ck;
476 dE(2) = -M_SQRT1_2 * (1.1 * c0 - 1.2 * c4) * ck;
485 double s0 = sin(
kappa * x(0) + 0.0 * M_PI);
486 double s4 = sin(
kappa * x(0) + 0.4 * M_PI);
487 double s9 = sin(
kappa * x(0) + 0.9 * M_PI);
489 f(0) = 2.2 * s0 + 1.2 * M_SQRT1_2 * s4;
491 M_SQRT1_2 * (1.1 * s0 + 1.3 * s9);
492 f(2) = 1.3 * (2.0 +
kappa *
kappa) * s9 + 1.2 * M_SQRT1_2 * s4;
496 double s0 = sin(
kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.0 * M_PI);
497 double s4 = sin(
kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.4 * M_PI);
498 double s9 = sin(
kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.9 * M_PI);
502 f(1) = 0.55 * (M_SQRT2 - kappa *
kappa) * s0 +
503 0.6 * (4.0 + kappa * kappa) * s4 +
505 f(2) = 0.6 * M_SQRT2 * s4 + 1.3 * (2.0 + kappa *
kappa) * s9;
509 double s0 = sin(kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.0 * M_PI);
510 double c0 = cos(kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.0 * M_PI);
511 double s4 = sin(kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.4 * M_PI);
512 double c4 = cos(kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.4 * M_PI);
513 double s9 = sin(kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.9 * M_PI);
514 double c9 = cos(kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.9 * M_PI);
515 double sk = sin(kappa * x(2));
516 double ck = cos(kappa * x(2));
518 f(0) = 0.55 * (4.0 + 3.0 * kappa *
kappa) * s0 * ck +
519 0.6 * (M_SQRT2 - kappa * kappa) * s4 * ck -
520 0.65 * M_SQRT2 * kappa * kappa * c9 * sk;
522 f(1) = 0.55 * (M_SQRT2 - kappa *
kappa) * s0 * ck +
523 0.6 * (4.0 + 3.0 * kappa * kappa) * s4 * ck +
524 0.65 * M_SQRT2 * s9 * ck -
525 0.65 * M_SQRT2 * kappa * kappa * c9 * sk;
527 f(2) = 0.6 * M_SQRT2 * s4 * ck -
528 M_SQRT2 * kappa * kappa * (0.55 * c0 + 0.6 * c4) * sk
529 + 1.3 * (2.0 + kappa * kappa) * s9 * ck;
A matrix coefficient that is constant in space and time.
int Size() const
Return the logical size of the array.
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get()...
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
The Auxiliary-space Maxwell Solver in hypre.
A coefficient that is constant across space and time.
Scalar coefficient defined as the inner product of two vector coefficients.
Arbitrary order 3D H(curl)-conforming Nedelec finite elements in 1D.
Integrator for (curl u, curl v) for Nedelec elements.
Vector coefficient that is constant in space and time.
Pointer to an Operator of a specified type.
HYPRE_BigInt GlobalTrueVSize() const
Data type dense matrix using column-major storage.
Arbitrary order 3D H(curl)-conforming Nedelec finite elements in 2D.
virtual void Save(std::ostream &out) const
Abstract parallel finite element space.
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
void SetPrintLevel(int print_lvl)
void SetOperator(const Operator &op)
Set/update the solver for the given operator.
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Vector coefficient defined as the Curl of a vector GridFunction.
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
void SetMaxIter(int max_iter)
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
A general vector function coefficient.
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
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...
void f_exact(const Vector &, Vector &)
Vector coefficient defined as a product of a matrix coefficient and a vector coefficient.
void SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
void E_exact(const Vector &, Vector &)
for VectorFiniteElements (Nedelec, Raviart-Thomas)
Arbitrary order H(curl)-conforming Nedelec finite elements.
Vector coefficient defined by a vector GridFunction.
void CurlE_exact(const Vector &, Vector &)
Arbitrary order H1-conforming (continuous) finite elements.
void Print(std::ostream &out=mfem::out) const override
Class for parallel grid function.
Wrapper for hypre's ParCSR matrix class.
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve Ax=b with hypre's PCG.
Class for parallel meshes.
void ParseCheck(std::ostream &out=mfem::out)
virtual double ComputeHCurlError(VectorCoefficient *exsol, VectorCoefficient *excurl, const IntegrationRule *irs[]=NULL) const
Returns the error measured H(curl)-norm for ND elements.
Arbitrary order "L2-conforming" discontinuous finite elements.
double f(const Vector &p)
double sigma(const Vector &x)