41 int main(
int argc,
char *argv[])
44 Mpi::Init(argc, argv);
45 int num_procs = Mpi::WorldSize();
46 int myid = Mpi::WorldRank();
50 const char *mesh_file =
"../data/inline-quad.mesh";
51 int ser_ref_levels = 2;
52 int par_ref_levels = 1;
55 bool visualization = 1;
58 args.
AddOption(&mesh_file,
"-m",
"--mesh",
60 args.
AddOption(&ser_ref_levels,
"-rs",
"--refine-serial",
61 "Number of times to refine the mesh uniformly in serial.");
62 args.
AddOption(&par_ref_levels,
"-rp",
"--refine-parallel",
63 "Number of times to refine the mesh uniformly in parallel.");
65 "Finite element order (polynomial degree) or -1 for"
66 " isoparametric space.");
68 "Number of desired eigenmodes.");
69 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
71 "Enable or disable GLVis visualization.");
77 Mesh *mesh =
new Mesh(mesh_file, 1, 1);
87 for (
int lev = 0; lev < ser_ref_levels; lev++)
96 ParMesh pmesh(MPI_COMM_WORLD, *mesh);
98 for (
int lev = 0; lev < par_ref_levels; lev++)
128 cout <<
"Number of H(Curl) unknowns: " << size_nd << endl;
129 cout <<
"Number of H(Div) unknowns: " << size_rt << endl;
146 epsilonMat(0,0) = 2.0; epsilonMat(1,1) = 2.0; epsilonMat(2,2) = 2.0;
147 epsilonMat(0,2) = 0.0; epsilonMat(2,0) = 0.0;
148 epsilonMat(0,1) = M_SQRT1_2; epsilonMat(1,0) = M_SQRT1_2;
149 epsilonMat(1,2) = M_SQRT1_2; epsilonMat(2,1) = M_SQRT1_2;
170 cout <<
"Computing eigenvalues shifted by " << shift << endl;
221 ostringstream mesh_name, mode_name, mode_deriv_name;
222 mesh_name <<
"mesh." << setfill(
'0') << setw(6) << myid;
224 ofstream mesh_ofs(mesh_name.str().c_str());
225 mesh_ofs.precision(8);
226 pmesh.
Print(mesh_ofs);
228 for (
int i=0; i<nev; i++)
234 mode_name <<
"mode_" << setfill(
'0') << setw(2) << i <<
"."
235 << setfill(
'0') << setw(6) << myid;
236 mode_deriv_name <<
"mode_deriv_" << setfill(
'0') << setw(2) << i <<
"."
237 << setfill(
'0') << setw(6) << myid;
239 ofstream mode_ofs(mode_name.str().c_str());
240 mode_ofs.precision(8);
244 ofstream mode_deriv_ofs(mode_deriv_name.str().c_str());
245 mode_deriv_ofs.precision(8);
246 dx.
Save(mode_deriv_ofs);
247 mode_deriv_name.str(
"");
263 mode_x_sock.precision(8);
264 mode_y_sock.precision(8);
265 mode_z_sock.precision(8);
266 mode_dy_sock.precision(8);
267 mode_dz_sock.precision(8);
269 Vector xVec(3); xVec = 0.0; xVec(0) = 1;
270 Vector yVec(3); yVec = 0.0; yVec(1) = 1;
271 Vector zVec(3); zVec = 0.0; zVec(2) = 1;
289 for (
int i=0; i<nev; i++)
293 cout <<
"Eigenmode " << i+1 <<
'/' << nev
294 <<
", Lambda = " << eigenvalues[i] - shift << endl;
314 double max_r = std::max(max_x, std::max(max_y, max_z));
317 x_cmd <<
" window_title 'Eigenmode " << i+1 <<
'/' << nev
318 <<
" X, Lambda = " << eigenvalues[i] - shift <<
"'"
319 <<
" valuerange -"<< max_r <<
' ' << max_r;
323 <<
" window_geometry 0 0 400 350";
326 mode_x_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n"
327 <<
"solution\n" << pmesh << xComp << flush
328 << x_cmd.str() << endl << flush;
330 MPI_Barrier(MPI_COMM_WORLD);
333 y_cmd <<
" window_title 'Eigenmode " << i+1 <<
'/' << nev
334 <<
" Y, Lambda = " << eigenvalues[i] - shift <<
"'"
335 <<
" valuerange -"<< max_r <<
' ' << max_r;
339 <<
" window_geometry 403 0 400 350";
342 mode_y_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n"
343 <<
"solution\n" << pmesh << yComp << flush
344 << y_cmd.str() << endl << flush;
346 MPI_Barrier(MPI_COMM_WORLD);
349 z_cmd <<
" window_title 'Eigenmode " << i+1 <<
'/' << nev
350 <<
" Z, Lambda = " << eigenvalues[i] - shift <<
"'"
351 <<
" valuerange -"<< max_r <<
' ' << max_r;
355 <<
" window_geometry 806 0 400 350";
358 mode_z_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n"
359 <<
"solution\n" << pmesh << zComp << flush
360 << z_cmd.str() << endl << flush;
362 MPI_Barrier(MPI_COMM_WORLD);
371 double min_d = max_r / (bbMax[0] - bbMin[0]);
375 max_r = std::max(std::max(max_y, max_z), min_d);
377 ostringstream dy_cmd;
378 dy_cmd <<
" window_title 'Curl Eigenmode "
380 <<
" Y, Lambda = " << eigenvalues[i] - shift <<
"'"
381 <<
"valuerange -"<< max_r <<
' ' << max_r;
385 <<
" window_geometry 403 375 400 350";
388 mode_dy_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n"
389 <<
"solution\n" << pmesh << dyComp << flush
390 << dy_cmd.str() << endl << flush;
392 MPI_Barrier(MPI_COMM_WORLD);
394 ostringstream dz_cmd;
395 dz_cmd <<
" window_title 'Curl Eigenmode "
397 <<
" Z, Lambda = " << eigenvalues[i] - shift <<
"'"
398 <<
"valuerange -"<< max_r <<
' ' << max_r;
402 <<
" window_geometry 806 375 400 350";
405 mode_dz_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n"
406 <<
"solution\n" << pmesh << dzComp << flush
407 << dz_cmd.str() << endl << flush;
409 MPI_Barrier(MPI_COMM_WORLD);
414 cout <<
"press (q)uit or (c)ontinue --> " << flush;
417 MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
427 mode_dy_sock.
close();
428 mode_dz_sock.
close();
436 mode_xy_sock.precision(8);
437 mode_z_sock.precision(8);
438 mode_dxy_sock.precision(8);
439 mode_dz_sock.precision(8);
442 xyMat(0,0) = 1.0; xyMat(1,1) = 1.0;
444 Vector zVec(3); zVec = 0.0; zVec(2) = 1;
463 for (
int i=0; i<nev; i++)
467 cout <<
"Eigenmode " << i+1 <<
'/' << nev
468 <<
", Lambda = " << eigenvalues[i] - shift << endl;
485 double max_r = std::max(max_v, max_s);
487 ostringstream xy_cmd;
488 xy_cmd <<
" window_title 'Eigenmode " << i+1 <<
'/' << nev
489 <<
" XY, Lambda = " << eigenvalues[i] - shift <<
"'"
490 <<
" valuerange 0.0 " << max_r;
493 xy_cmd <<
" keys aavvv"
494 <<
" window_geometry 0 0 400 350";
497 mode_xy_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n"
498 <<
"solution\n" << pmesh << xyComp << flush
499 << xy_cmd.str() << endl << flush;
501 MPI_Barrier(MPI_COMM_WORLD);
504 z_cmd <<
" window_title 'Eigenmode " << i+1 <<
'/' << nev
505 <<
" Z, Lambda = " << eigenvalues[i] - shift <<
"'"
506 <<
" valuerange -"<< max_r <<
' ' << max_r;
510 <<
" window_geometry 403 0 400 350";
513 mode_z_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n"
514 <<
"solution\n" << pmesh << zComp << flush
515 << z_cmd.str() << endl << flush;
517 MPI_Barrier(MPI_COMM_WORLD);
526 double min_d = max_r / std::min(bbMax[0] - bbMin[0],
527 bbMax[1] - bbMin[1]);
531 max_r = std::max(std::max(max_v, max_s), min_d);
533 ostringstream dxy_cmd;
534 dxy_cmd <<
" window_title 'Curl Eigenmode "
536 <<
" XY, Lambda = " << eigenvalues[i] - shift <<
"'"
537 <<
" valuerange 0.0 " << max_r <<
'\n';
540 dxy_cmd <<
" keys aavvv "
541 <<
" window_geometry 0 375 400 350";
545 mode_dxy_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n"
546 <<
"solution\n" << pmesh << dxyComp << flush
547 << dxy_cmd.str() << endl << flush;
549 MPI_Barrier(MPI_COMM_WORLD);
551 ostringstream dz_cmd;
552 dz_cmd <<
" window_title 'Curl Eigenmode "
554 <<
" Z, Lambda = " << eigenvalues[i] - shift <<
"'"
555 <<
" valuerange -" << max_r <<
' ' << max_r;
559 <<
" window_geometry 403 375 400 350";
562 mode_dz_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n"
563 <<
"solution\n" << pmesh << dzComp << flush
564 << dz_cmd.str() << endl << flush;
566 MPI_Barrier(MPI_COMM_WORLD);
571 cout <<
"press (q)uit or (c)ontinue --> " << flush;
574 MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
581 mode_xy_sock.
close();
583 mode_dxy_sock.
close();
584 mode_dz_sock.
close();
590 mode_sock.precision(8);
591 mode_deriv_sock.precision(8);
593 for (
int i=0; i<nev; i++)
597 cout <<
"Eigenmode " << i+1 <<
'/' << nev
598 <<
", Lambda = " << eigenvalues[i] - shift << endl;
605 mode_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n"
606 <<
"solution\n" << pmesh << x << flush
607 <<
"window_title 'Eigenmode " << i+1 <<
'/' << nev
608 <<
", Lambda = " << eigenvalues[i] - shift
611 MPI_Barrier(MPI_COMM_WORLD);
613 mode_deriv_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n"
614 <<
"solution\n" << pmesh << dx << flush
615 <<
"window_geometry 0 375 400 350 "
616 <<
"window_title 'Curl Eigenmode "
618 <<
", Lambda = " << eigenvalues[i] - shift
621 MPI_Barrier(MPI_COMM_WORLD);
626 cout <<
"press (q)uit or (c)ontinue --> " << flush;
629 MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
654 Vector zeroVec(vdim); zeroVec = 0.0;
A matrix coefficient that is constant in space and time.
int Size() const
Return the logical size of the array.
The Auxiliary-space Maxwell Solver in hypre.
A coefficient that is constant across space and time.
void SetMassMatrix(const HypreParMatrix &M)
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.
double GetVectorMax(int vdim, const ParGridFunction &x)
Vector coefficient that is constant in space and time.
void GetBoundingBox(Vector &min, Vector &max, int ref=2)
Returns the minimum and maximum corners of the mesh bounding box.
HYPRE_BigInt GlobalTrueVSize() const
void SetPreconditioner(HypreSolver &precond)
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 AddDomainInterpolator(DiscreteInterpolator *di)
Adds a domain interpolator. Assumes ownership of di.
int close()
Close the socketstream.
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
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 SetPrintLevel(int print_lvl)
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 SetSize(int nsize)
Change the logical size of the array, keep existing entries.
void Solve()
Solve the eigenproblem.
virtual double ComputeMaxError(Coefficient *exsol[], const IntegrationRule *irs[]=NULL) const
void SetNumModes(int num_eigs)
const HypreParVector & GetEigenvector(unsigned int i) const
Extract a single eigenvector.
Vector coefficient defined as a product of a matrix coefficient and a vector coefficient.
void SetMaxIter(int max_iter)
void GetEigenvalues(Array< double > &eigenvalues) const
Collect the converged eigenvalues.
void SetPrintLevel(int logging)
Arbitrary order H(curl)-conforming Nedelec finite elements.
Vector coefficient defined by a vector GridFunction.
void SetOperator(const HypreParMatrix &A)
Arbitrary order 3D H(div)-conforming Raviart-Thomas finite elements in 2D.
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 Assemble(int skip_zeros=1)
Construct the internal matrix representation of the discrete linear operator.
Class for parallel meshes.
void SetSingularProblem()
Set this option when solving a curl-curl problem with zero mass term.
Arbitrary order 3D H(div)-conforming Raviart-Thomas finite elements in 1D.
void ParseCheck(std::ostream &out=mfem::out)
double GetScalarMax(const ParGridFunction &x)
Arbitrary order "L2-conforming" discontinuous finite elements.