35 #define jn(n, x) _jn(n, x)
36 #define yn(n, x) _yn(n, x)
77 Array<int> * GetMarkedPMLElements() {
return &elems;}
80 void SetAttributes(
ParMesh *pmesh);
83 void StretchFunction(
const Vector &x, vector<complex<double>> &dxs);
90 CartesianPML * pml =
nullptr;
91 void (*Function)(
const Vector &, CartesianPML * ,
Vector &);
93 PMLDiagMatrixCoefficient(
int dim,
void(*F)(
const Vector &, CartesianPML *,
99 using VectorCoefficient::Eval;
108 (*Function)(transip, pml, K);
151 int main(
int argc,
char *argv[])
155 MPI_Init(&argc, &argv);
156 MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
157 MPI_Comm_rank(MPI_COMM_WORLD, &myid);
160 const char *mesh_file =
nullptr;
163 int par_ref_levels = 2;
166 bool herm_conv =
true;
167 bool visualization = 1;
170 args.
AddOption(&mesh_file,
"-m",
"--mesh",
171 "Mesh file to use.");
173 "Finite element order (polynomial degree).");
174 args.
AddOption(&iprob,
"-prob",
"--problem",
"Problem case"
175 " 0: beam, 1: disc, 2: lshape, 3: fichera, 4: General");
176 args.
AddOption(&ref_levels,
"-rs",
"--refinements-serial",
177 "Number of serial refinements");
178 args.
AddOption(&par_ref_levels,
"-rp",
"--refinements-parallel",
179 "Number of parallel refinements");
181 "Permeability of free space (or 1/(spring constant)).");
183 "Permittivity of free space (or mass constant).");
184 args.
AddOption(&freq,
"-f",
"--frequency",
185 "Frequency (in Hz).");
186 args.
AddOption(&herm_conv,
"-herm",
"--hermitian",
"-no-herm",
187 "--no-hermitian",
"Use convention for Hermitian operators.");
188 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
189 "--no-visualization",
190 "Enable or disable GLVis visualization.");
193 if (iprob > 4) { iprob = 4; }
203 mesh_file =
"../data/beam-hex.mesh";
206 mesh_file =
"../data/square-disc.mesh";
209 mesh_file =
"../data/l-shape.mesh";
212 mesh_file =
"../data/fichera.mesh";
216 mesh_file =
"../data/inline-quad.mesh";
235 Mesh * mesh =
new Mesh(mesh_file, 1, 1);
266 CartesianPML * pml =
new CartesianPML(mesh,length);
271 for (
int l = 0; l < ref_levels; l++)
280 for (
int l = 0; l < par_ref_levels; l++)
290 pml->SetAttributes(pmesh);
299 cout <<
"Number of finite element unknowns: " << size << endl;
314 for (
int j = 0; j < pmesh->
GetNBE(); j++)
324 if (center[0] == 1.0 || center[0] == 0.5 || center[1] == 0.5)
330 if (center[0] == -1.0 || center[0] == 0.0 ||
331 center[1] == 0.0 || center[2] == 0.0)
346 herm_conv ? ComplexOperator::HERMITIAN : ComplexOperator::BLOCK_SYMMETRIC;
356 b.Vector::operator=(0.0);
382 attr = 0; attr[0] = 1;
400 int cdim = (dim == 2) ? 1 : dim;
432 #ifdef MFEM_USE_SUPERLU
488 (conv == ComplexOperator::HERMITIAN) ? -1.0 : 1.0);
513 int order_quad = max(2, 2 * order + 1);
515 for (
int i = 0; i < Geometry::NumGeom; ++i)
521 pml->GetMarkedPMLElements());
523 pml->GetMarkedPMLElements());
527 double norm_E_Re, norm_E_Im;
529 pml->GetMarkedPMLElements());
531 pml->GetMarkedPMLElements());
535 cout <<
"\n Relative Error (Re part): || E_h - E || / ||E|| = "
536 << L2Error_Re / norm_E_Re
537 <<
"\n Relative Error (Im part): || E_h - E || / ||E|| = "
538 << L2Error_Im / norm_E_Im
539 <<
"\n Total Error: "
540 << sqrt(L2Error_Re*L2Error_Re + L2Error_Im*L2Error_Im) <<
"\n\n";
547 ostringstream mesh_name, sol_r_name, sol_i_name;
548 mesh_name <<
"mesh." << setfill(
'0') << setw(6) << myid;
549 sol_r_name <<
"ex25p-sol_r." << setfill(
'0') << setw(6) << myid;
550 sol_i_name <<
"ex25p-sol_i." << setfill(
'0') << setw(6) << myid;
552 ofstream mesh_ofs(mesh_name.str().c_str());
553 mesh_ofs.precision(8);
554 pmesh->
Print(mesh_ofs);
556 ofstream sol_r_ofs(sol_r_name.str().c_str());
557 ofstream sol_i_ofs(sol_i_name.str().c_str());
558 sol_r_ofs.precision(8);
559 sol_i_ofs.precision(8);
569 keys = (dim == 3) ?
"keys macF\n" : keys =
"keys amrRljcUUuu\n";
570 if (
prob ==
beam && dim == 3) {keys =
"keys macFFiYYYYYYYYYYYYYYYYYY\n";}
571 if (
prob ==
beam && dim == 2) {keys =
"keys amrRljcUUuuu\n"; }
578 sol_sock_re.precision(8);
579 sol_sock_re <<
"parallel " << num_procs <<
" " << myid <<
"\n"
580 <<
"solution\n" << *pmesh << x.
real() << keys
581 <<
"window_title 'Solution real part'" << flush;
582 MPI_Barrier(MPI_COMM_WORLD);
587 sol_sock_im.precision(8);
588 sol_sock_im <<
"parallel " << num_procs <<
" " << myid <<
"\n"
589 <<
"solution\n" << *pmesh << x.
imag() << keys
590 <<
"window_title 'Solution imag part'" << flush;
591 MPI_Barrier(MPI_COMM_WORLD);
599 sol_sock.precision(8);
600 sol_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n"
601 <<
"solution\n" << *pmesh << x_t << keys <<
"autoscale off\n"
602 <<
"window_title 'Harmonic Solution (t = 0.0 T)'"
603 <<
"pause\n" << flush;
607 cout <<
"GLVis visualization paused."
608 <<
" Press space (in the GLVis window) to resume it.\n";
615 double t = (double)(i % num_frames) / num_frames;
617 oss <<
"Harmonic Solution (t = " << t <<
" T)";
619 add(cos(2.0*M_PI*t), x.
real(), sin(2.0*M_PI*t), x.
imag(), x_t);
620 sol_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n";
621 sol_sock <<
"solution\n" << *pmesh << x_t
622 <<
"window_title '" << oss.str() <<
"'" << flush;
641 for (
int i = 0; i <
dim; ++i)
644 r += pow(x[i] - center[i], 2.);
647 double coeff = pow(n, 2) / M_PI;
648 double alpha = -pow(n, 2) * r;
650 f[0] = coeff * exp(alpha);
656 for (
int i = 0; i <
dim; ++i)
661 complex<double> zi = complex<double>(0., 1.);
677 double x0 = x(0) + shift(0);
678 double x1 = x(1) + shift(1);
679 double r = sqrt(x0 * x0 + x1 * x1);
683 complex<double> Ho, Ho_r, Ho_rr;
684 Ho = jn(0, beta) + zi * yn(0, beta);
685 Ho_r = -k * (jn(1, beta) + zi * yn(1, beta));
686 Ho_rr = -k * k * (1.0 / beta *
687 (jn(1, beta) + zi * yn(1, beta)) -
688 (jn(2, beta) + zi * yn(2, beta)));
693 double r_xy = -(r_x / r) * r_y;
694 double r_xx = (1.0 / r) * (1.0 - r_x * r_x);
696 complex<double> val, val_xx, val_xy;
697 val = 0.25 * zi * Ho;
698 val_xx = 0.25 * zi * (r_xx * Ho_r + r_x * r_x * Ho_rr);
699 val_xy = 0.25 * zi * (r_xy * Ho_r + r_x * r_y * Ho_rr);
700 E[0] = zi / k * (k * k * val + val_xx);
701 E[1] = zi / k * val_xy;
705 double x0 = x(0) + shift(0);
706 double x1 = x(1) + shift(1);
707 double x2 = x(2) + shift(2);
708 double r = sqrt(x0 * x0 + x1 * x1 + x2 * x2);
713 double r_xx = (1.0 / r) * (1.0 - r_x * r_x);
714 double r_yx = -(r_y / r) * r_x;
715 double r_zx = -(r_z / r) * r_x;
717 complex<double> val, val_r, val_rr;
718 val = exp(zi * k * r) / r;
719 val_r = val / r * (zi * k * r - 1.0);
720 val_rr = val / (r * r) * (-k * k * r * r
721 - 2.0 * zi * k * r + 2.0);
723 complex<double> val_xx, val_yx, val_zx;
724 val_xx = val_rr * r_x * r_x + val_r * r_xx;
725 val_yx = val_rr * r_x * r_y + val_r * r_yx;
726 val_zx = val_rr * r_x * r_z + val_r * r_zx;
728 complex<double>
alpha = zi * k / 4.0 / M_PI / k / k;
729 E[0] = alpha * (k * k * val + val_xx);
730 E[1] = alpha * val_yx;
731 E[2] = alpha * val_zx;
740 double k10 = sqrt(k * k - M_PI * M_PI);
741 E[1] = -zi * k / M_PI * sin(M_PI*x(2))*exp(zi * k10 * x(0));
745 E[1] = -zi * k / M_PI * exp(zi * k * x(0));
756 vector<complex<double>> Eval(E.
Size());
758 for (
int i = 0; i <
dim; ++i)
760 E[i] = Eval[i].real();
766 vector<complex<double>> Eval(E.
Size());
768 for (
int i = 0; i <
dim; ++i)
770 E[i] = Eval[i].imag();
779 for (
int i = 0; i <
dim; ++i)
791 vector<complex<double>> Eval(E.
Size());
793 for (
int i = 0; i <
dim; ++i)
795 E[i] = Eval[i].real();
806 for (
int i = 0; i <
dim; ++i)
818 vector<complex<double>> Eval(E.
Size());
820 for (
int i = 0; i <
dim; ++i)
822 E[i] = Eval[i].imag();
829 vector<complex<double>> dxs(dim);
830 complex<double> det(1.0, 0.0);
831 pml->StretchFunction(x, dxs);
833 for (
int i = 0; i <
dim; ++i)
838 for (
int i = 0; i <
dim; ++i)
840 D(i) = (det / pow(dxs[i], 2)).real();
846 vector<complex<double>> dxs(dim);
847 complex<double> det = 1.0;
848 pml->StretchFunction(x, dxs);
850 for (
int i = 0; i <
dim; ++i)
855 for (
int i = 0; i <
dim; ++i)
857 D(i) = (det / pow(dxs[i], 2)).imag();
863 vector<complex<double>> dxs(dim);
864 complex<double> det = 1.0;
865 pml->StretchFunction(x, dxs);
867 for (
int i = 0; i <
dim; ++i)
872 for (
int i = 0; i <
dim; ++i)
874 D(i) = abs(det / pow(dxs[i], 2));
880 vector<complex<double>> dxs(dim);
881 complex<double> det(1.0, 0.0);
882 pml->StretchFunction(x, dxs);
884 for (
int i = 0; i <
dim; ++i)
892 D = (1.0 / det).real();
896 for (
int i = 0; i <
dim; ++i)
898 D(i) = (pow(dxs[i], 2) / det).real();
905 vector<complex<double>> dxs(dim);
906 complex<double> det = 1.0;
907 pml->StretchFunction(x, dxs);
909 for (
int i = 0; i <
dim; ++i)
916 D = (1.0 / det).imag();
920 for (
int i = 0; i <
dim; ++i)
922 D(i) = (pow(dxs[i], 2) / det).imag();
929 vector<complex<double>> dxs(dim);
930 complex<double> det = 1.0;
931 pml->StretchFunction(x, dxs);
933 for (
int i = 0; i <
dim; ++i)
944 for (
int i = 0; i <
dim; ++i)
946 D(i) = abs(pow(dxs[i], 2) / det);
952 : mesh(mesh_), length(length_)
954 dim = mesh->Dimension();
958 void CartesianPML::SetBoundaries()
960 comp_dom_bdr.SetSize(dim, 2);
961 dom_bdr.SetSize(dim, 2);
963 mesh->GetBoundingBox(pmin, pmax);
964 for (
int i = 0; i <
dim; i++)
966 dom_bdr(i, 0) = pmin(i);
967 dom_bdr(i, 1) = pmax(i);
968 comp_dom_bdr(i, 0) = dom_bdr(i, 0) + length(i, 0);
969 comp_dom_bdr(i, 1) = dom_bdr(i, 1) - length(i, 1);
973 void CartesianPML::SetAttributes(
ParMesh *pmesh)
976 MPI_Comm_rank(MPI_COMM_WORLD,&myid);
977 int nrelem = pmesh->
GetNE();
980 elems.SetSize(nrelem);
983 for (
int i = 0; i < nrelem; ++i)
993 int nrvert = vertices.
Size();
996 for (
int iv = 0; iv < nrvert; ++iv)
998 int vert_idx = vertices[iv];
999 double *coords = pmesh->
GetVertex(vert_idx);
1000 for (
int comp = 0; comp <
dim; ++comp)
1002 if (coords[comp] > comp_dom_bdr(comp, 1) ||
1003 coords[comp] < comp_dom_bdr(comp, 0))
1019 void CartesianPML::StretchFunction(
const Vector &x,
1020 vector<complex<double>> &dxs)
1022 complex<double> zi = complex<double>(0., 1.);
1030 for (
int i = 0; i <
dim; ++i)
1035 coeff = n * c / k / pow(length(i, 1), n);
1036 dxs[i] = 1.0 + zi * coeff *
1041 coeff = n * c / k / pow(length(i, 0), n);
1042 dxs[i] = 1.0 + zi * coeff *
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()...
Class for an integration rule - an Array of IntegrationPoint.
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
int GetBdrAttribute(int i) const
Return the attribute of boundary element i.
The Auxiliary-space Maxwell Solver in hypre.
const double * GetVertex(int i) const
Return pointer to vertex i's coordinates.
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Base class for vector Coefficients that optionally depend on time and space.
A coefficient that is constant across space and time.
virtual void GetVertices(Array< int > &v) const =0
Returns element's vertices.
void E_exact_Re(const Vector &x, Vector &E)
virtual void ReorientTetMesh()
See the remarks for the serial version in mesh.hpp.
void SetSize(int s)
Resize the vector to size s.
int GetNBE() const
Returns number of boundary elements.
Integrator for (curl u, curl v) for Nedelec elements.
Pointer to an Operator of a specified type.
int Size() const
Returns the size of the vector.
void maxwell_solution(const Vector &x, vector< complex< double >> &Eval)
int GetNE() const
Returns number of elements.
void detJ_JT_J_inv_abs(const Vector &x, CartesianPML *pml, Vector &D)
virtual void Save(std::ostream &out) const
Abstract parallel finite element space.
void detJ_inv_JT_J_abs(const Vector &x, CartesianPML *pml, Vector &D)
void E_bdr_data_Re(const Vector &x, Vector &E)
virtual void SetAttributes()
const IntegrationPoint & GetCenter(int GeomType)
Return the center of the given Geometry::Type, GeomType.
int main(int argc, char *argv[])
void add(const Vector &v1, const Vector &v2, Vector &v)
void SetSymmetricPattern(bool sym)
virtual int GetTrueVSize() const
Return the number of local vector true dofs.
ElementTransformation * GetBdrElementTransformation(int i)
Returns the transformation defining the i-th boundary element.
void SetOperator(const Operator &op)
Set/update the solver for the given operator.
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Derived vector coefficient that has the value of the parent vector where it is active and is zero oth...
void source(const Vector &x, Vector &f)
void E_exact_Im(const Vector &x, Vector &E)
void SetPrintLevel(int print_lvl)
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Array2D< double > domain_bdr
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
A class to handle Block diagonal preconditioners in a matrix-free implementation. ...
Derived coefficient that takes the value of the parent coefficient for the active attributes and is z...
void E_bdr_data_Im(const Vector &x, Vector &E)
void SetColumnPermutation(superlu::ColPerm col_perm)
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
void detJ_inv_JT_J_Im(const Vector &x, CartesianPML *pml, Vector &D)
void SetMaxIter(int max_it)
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
HYPRE_Int GlobalTrueVSize() const
void SetKDim(int dim)
Set the number of iteration to perform between restarts, default is 50.
virtual void Print(std::ostream &out=mfem::out) const
void SetPrintStatistics(bool print_stat)
const Element * GetElement(int i) const
void PrintUsage(std::ostream &out) const
Print the usage message.
A general vector function coefficient.
void SetAbsTol(double atol)
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Array2D< double > comp_domain_bdr
void SetRelTol(double rtol)
Scaled Operator B: x -> a A(x).
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Vector coefficient defined as a product of scalar and vector coefficients.
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 PartialSum()
Fill the entries of the array with the cumulative sum of the entries.
void detJ_JT_J_inv_Re(const Vector &x, CartesianPML *pml, Vector &D)
virtual void ProjectBdrCoefficientTangent(VectorCoefficient &real_coeff, VectorCoefficient &imag_coeff, Array< int > &attr)
Class for integration point with weight.
virtual double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL) const
void detJ_JT_J_inv_Im(const Vector &x, CartesianPML *pml, Vector &D)
Specialization of the ComplexOperator built from a pair of HypreParMatrices.
void PrintOptions(std::ostream &out) const
Print the options.
void SetAttribute(const int attr)
Set element's attribute.
Geometry::Type GetBdrElementBaseGeometry(int i) const
for VectorFiniteElements (Nedelec, Raviart-Thomas)
Arbitrary order H(curl)-conforming Nedelec finite elements.
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Class for parallel grid function.
void detJ_inv_JT_J_Re(const Vector &x, CartesianPML *pml, Vector &D)
Wrapper for hypre's ParCSR matrix class.
Class for parallel meshes.
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Abstract data type element.
Array< int > attributes
A list of all unique element attributes used by the Mesh.
void SetDiagonalBlock(int iblock, Operator *op)
Add a square block op in the block-entry (iblock, iblock).
double f(const Vector &p)
bool Good() const
Return true if the command line options were parsed successfully.