39 #define jn(n, x) _jn(n, x)
40 #define yn(n, x) _yn(n, x)
81 Array<int> * GetMarkedPMLElements() {
return &elems;}
84 void SetAttributes(
ParMesh *pmesh);
87 void StretchFunction(
const Vector &x, vector<complex<double>> &dxs);
94 CartesianPML * pml =
nullptr;
95 void (*Function)(
const Vector &, CartesianPML * ,
Vector &);
97 PMLDiagMatrixCoefficient(
int dim,
void(*F)(
const Vector &, CartesianPML *,
103 using VectorCoefficient::Eval;
112 (*Function)(transip, pml, K);
155 int main(
int argc,
char *argv[])
159 MPI_Init(&argc, &argv);
160 MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
161 MPI_Comm_rank(MPI_COMM_WORLD, &myid);
164 const char *mesh_file =
nullptr;
167 int par_ref_levels = 2;
170 bool herm_conv =
true;
171 bool slu_solver =
false;
172 bool mumps_solver =
false;
173 bool visualization = 1;
175 const char *device_config =
"cpu";
178 args.
AddOption(&mesh_file,
"-m",
"--mesh",
179 "Mesh file to use.");
181 "Finite element order (polynomial degree).");
182 args.
AddOption(&iprob,
"-prob",
"--problem",
"Problem case"
183 " 0: beam, 1: disc, 2: lshape, 3: fichera, 4: General");
184 args.
AddOption(&ref_levels,
"-rs",
"--refinements-serial",
185 "Number of serial refinements");
186 args.
AddOption(&par_ref_levels,
"-rp",
"--refinements-parallel",
187 "Number of parallel refinements");
189 "Permeability of free space (or 1/(spring constant)).");
191 "Permittivity of free space (or mass constant).");
192 args.
AddOption(&freq,
"-f",
"--frequency",
193 "Frequency (in Hz).");
194 args.
AddOption(&herm_conv,
"-herm",
"--hermitian",
"-no-herm",
195 "--no-hermitian",
"Use convention for Hermitian operators.");
196 #ifdef MFEM_USE_SUPERLU
197 args.
AddOption(&slu_solver,
"-slu",
"--superlu",
"-no-slu",
198 "--no-superlu",
"Use the SuperLU Solver.");
200 #ifdef MFEM_USE_MUMPS
201 args.
AddOption(&mumps_solver,
"-mumps",
"--mumps-solver",
"-no-mumps",
202 "--no-mumps-solver",
"Use the MUMPS Solver.");
204 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
205 "--no-visualization",
206 "Enable or disable GLVis visualization.");
207 args.
AddOption(&pa,
"-pa",
"--partial-assembly",
"-no-pa",
208 "--no-partial-assembly",
"Enable Partial Assembly.");
209 args.
AddOption(&device_config,
"-d",
"--device",
210 "Device configuration string, see Device::Configure().");
212 if (slu_solver && mumps_solver)
215 cout <<
"WARNING: Both SuperLU and MUMPS have been selected,"
216 <<
" please choose either one." << endl
217 <<
" Defaulting to SuperLU." << endl;
218 mumps_solver =
false;
221 if (iprob > 4) { iprob = 4; }
226 Device device(device_config);
227 if (myid == 0) { device.
Print(); }
236 mesh_file =
"../data/beam-hex.mesh";
239 mesh_file =
"../data/square-disc.mesh";
242 mesh_file =
"../data/l-shape.mesh";
245 mesh_file =
"../data/fichera.mesh";
249 mesh_file =
"../data/inline-quad.mesh";
268 Mesh * mesh =
new Mesh(mesh_file, 1, 1);
299 CartesianPML * pml =
new CartesianPML(mesh,length);
304 for (
int l = 0; l < ref_levels; l++)
313 for (
int l = 0; l < par_ref_levels; l++)
323 pml->SetAttributes(pmesh);
332 cout <<
"Number of finite element unknowns: " << size << endl;
347 for (
int j = 0; j < pmesh->
GetNBE(); j++)
357 if (center[0] == 1.0 || center[0] == 0.5 || center[1] == 0.5)
363 if (center[0] == -1.0 || center[0] == 0.0 ||
364 center[1] == 0.0 || center[2] == 0.0)
379 herm_conv ? ComplexOperator::HERMITIAN : ComplexOperator::BLOCK_SYMMETRIC;
389 b.Vector::operator=(0.0);
415 attr = 0; attr[0] = 1;
433 int cdim = (dim == 2) ? 1 : dim;
466 #ifdef MFEM_USE_SUPERLU
467 if (!pa && slu_solver)
481 #ifdef MFEM_USE_MUMPS
482 if (!pa && mumps_solver)
500 if (pa || (!slu_solver && !mumps_solver))
533 int s = (conv == ComplexOperator::HERMITIAN) ? -1.0 : 1.0;
578 int order_quad = max(2, 2 * order + 1);
580 for (
int i = 0; i < Geometry::NumGeom; ++i)
586 pml->GetMarkedPMLElements());
588 pml->GetMarkedPMLElements());
592 double norm_E_Re, norm_E_Im;
594 pml->GetMarkedPMLElements());
596 pml->GetMarkedPMLElements());
600 cout <<
"\n Relative Error (Re part): || E_h - E || / ||E|| = "
601 << L2Error_Re / norm_E_Re
602 <<
"\n Relative Error (Im part): || E_h - E || / ||E|| = "
603 << L2Error_Im / norm_E_Im
604 <<
"\n Total Error: "
605 << sqrt(L2Error_Re*L2Error_Re + L2Error_Im*L2Error_Im) <<
"\n\n";
612 ostringstream mesh_name, sol_r_name, sol_i_name;
613 mesh_name <<
"mesh." << setfill(
'0') << setw(6) << myid;
614 sol_r_name <<
"ex25p-sol_r." << setfill(
'0') << setw(6) << myid;
615 sol_i_name <<
"ex25p-sol_i." << setfill(
'0') << setw(6) << myid;
617 ofstream mesh_ofs(mesh_name.str().c_str());
618 mesh_ofs.precision(8);
619 pmesh->
Print(mesh_ofs);
621 ofstream sol_r_ofs(sol_r_name.str().c_str());
622 ofstream sol_i_ofs(sol_i_name.str().c_str());
623 sol_r_ofs.precision(8);
624 sol_i_ofs.precision(8);
634 keys = (dim == 3) ?
"keys macF\n" : keys =
"keys amrRljcUUuu\n";
635 if (
prob ==
beam && dim == 3) {keys =
"keys macFFiYYYYYYYYYYYYYYYYYY\n";}
636 if (
prob ==
beam && dim == 2) {keys =
"keys amrRljcUUuuu\n"; }
643 sol_sock_re.precision(8);
644 sol_sock_re <<
"parallel " << num_procs <<
" " << myid <<
"\n"
645 <<
"solution\n" << *pmesh << x.
real() << keys
646 <<
"window_title 'Solution real part'" << flush;
647 MPI_Barrier(MPI_COMM_WORLD);
652 sol_sock_im.precision(8);
653 sol_sock_im <<
"parallel " << num_procs <<
" " << myid <<
"\n"
654 <<
"solution\n" << *pmesh << x.
imag() << keys
655 <<
"window_title 'Solution imag part'" << flush;
656 MPI_Barrier(MPI_COMM_WORLD);
664 sol_sock.precision(8);
665 sol_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n"
666 <<
"solution\n" << *pmesh << x_t << keys <<
"autoscale off\n"
667 <<
"window_title 'Harmonic Solution (t = 0.0 T)'"
668 <<
"pause\n" << flush;
672 cout <<
"GLVis visualization paused."
673 <<
" Press space (in the GLVis window) to resume it.\n";
680 double t = (double)(i % num_frames) / num_frames;
682 oss <<
"Harmonic Solution (t = " << t <<
" T)";
684 add(cos(2.0*M_PI*t), x.
real(), sin(2.0*M_PI*t), x.
imag(), x_t);
685 sol_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n";
686 sol_sock <<
"solution\n" << *pmesh << x_t
687 <<
"window_title '" << oss.str() <<
"'" << flush;
706 for (
int i = 0; i <
dim; ++i)
709 r += pow(x[i] - center[i], 2.);
712 double coeff = pow(n, 2) / M_PI;
713 double alpha = -pow(n, 2) * r;
715 f[0] = coeff * exp(alpha);
721 for (
int i = 0; i <
dim; ++i)
726 complex<double> zi = complex<double>(0., 1.);
742 double x0 = x(0) + shift(0);
743 double x1 = x(1) + shift(1);
744 double r = sqrt(x0 * x0 + x1 * x1);
748 complex<double> Ho, Ho_r, Ho_rr;
749 Ho = jn(0, beta) + zi * yn(0, beta);
750 Ho_r = -k * (jn(1, beta) + zi * yn(1, beta));
751 Ho_rr = -k * k * (1.0 / beta *
752 (jn(1, beta) + zi * yn(1, beta)) -
753 (jn(2, beta) + zi * yn(2, beta)));
758 double r_xy = -(r_x / r) * r_y;
759 double r_xx = (1.0 / r) * (1.0 - r_x * r_x);
761 complex<double> val, val_xx, val_xy;
762 val = 0.25 * zi * Ho;
763 val_xx = 0.25 * zi * (r_xx * Ho_r + r_x * r_x * Ho_rr);
764 val_xy = 0.25 * zi * (r_xy * Ho_r + r_x * r_y * Ho_rr);
765 E[0] = zi / k * (k * k * val + val_xx);
766 E[1] = zi / k * val_xy;
770 double x0 = x(0) + shift(0);
771 double x1 = x(1) + shift(1);
772 double x2 = x(2) + shift(2);
773 double r = sqrt(x0 * x0 + x1 * x1 + x2 * x2);
778 double r_xx = (1.0 / r) * (1.0 - r_x * r_x);
779 double r_yx = -(r_y / r) * r_x;
780 double r_zx = -(r_z / r) * r_x;
782 complex<double> val, val_r, val_rr;
783 val = exp(zi * k * r) / r;
784 val_r = val / r * (zi * k * r - 1.0);
785 val_rr = val / (r * r) * (-k * k * r * r
786 - 2.0 * zi * k * r + 2.0);
788 complex<double> val_xx, val_yx, val_zx;
789 val_xx = val_rr * r_x * r_x + val_r * r_xx;
790 val_yx = val_rr * r_x * r_y + val_r * r_yx;
791 val_zx = val_rr * r_x * r_z + val_r * r_zx;
793 complex<double>
alpha = zi * k / 4.0 / M_PI / k / k;
794 E[0] = alpha * (k * k * val + val_xx);
795 E[1] = alpha * val_yx;
796 E[2] = alpha * val_zx;
805 double k10 = sqrt(k * k - M_PI * M_PI);
806 E[1] = -zi * k / M_PI * sin(M_PI*x(2))*exp(zi * k10 * x(0));
810 E[1] = -zi * k / M_PI * exp(zi * k * x(0));
821 vector<complex<double>> Eval(E.
Size());
823 for (
int i = 0; i <
dim; ++i)
825 E[i] = Eval[i].real();
831 vector<complex<double>> Eval(E.
Size());
833 for (
int i = 0; i <
dim; ++i)
835 E[i] = Eval[i].imag();
844 for (
int i = 0; i <
dim; ++i)
856 vector<complex<double>> Eval(E.
Size());
858 for (
int i = 0; i <
dim; ++i)
860 E[i] = Eval[i].real();
871 for (
int i = 0; i <
dim; ++i)
883 vector<complex<double>> Eval(E.
Size());
885 for (
int i = 0; i <
dim; ++i)
887 E[i] = Eval[i].imag();
894 vector<complex<double>> dxs(dim);
895 complex<double> det(1.0, 0.0);
896 pml->StretchFunction(x, dxs);
898 for (
int i = 0; i <
dim; ++i)
903 for (
int i = 0; i <
dim; ++i)
905 D(i) = (det / pow(dxs[i], 2)).real();
911 vector<complex<double>> dxs(dim);
912 complex<double> det = 1.0;
913 pml->StretchFunction(x, dxs);
915 for (
int i = 0; i <
dim; ++i)
920 for (
int i = 0; i <
dim; ++i)
922 D(i) = (det / pow(dxs[i], 2)).imag();
928 vector<complex<double>> dxs(dim);
929 complex<double> det = 1.0;
930 pml->StretchFunction(x, dxs);
932 for (
int i = 0; i <
dim; ++i)
937 for (
int i = 0; i <
dim; ++i)
939 D(i) = abs(det / pow(dxs[i], 2));
945 vector<complex<double>> dxs(dim);
946 complex<double> det(1.0, 0.0);
947 pml->StretchFunction(x, dxs);
949 for (
int i = 0; i <
dim; ++i)
957 D = (1.0 / det).real();
961 for (
int i = 0; i <
dim; ++i)
963 D(i) = (pow(dxs[i], 2) / det).real();
970 vector<complex<double>> dxs(dim);
971 complex<double> det = 1.0;
972 pml->StretchFunction(x, dxs);
974 for (
int i = 0; i <
dim; ++i)
981 D = (1.0 / det).imag();
985 for (
int i = 0; i <
dim; ++i)
987 D(i) = (pow(dxs[i], 2) / det).imag();
994 vector<complex<double>> dxs(dim);
995 complex<double> det = 1.0;
996 pml->StretchFunction(x, dxs);
998 for (
int i = 0; i <
dim; ++i)
1009 for (
int i = 0; i <
dim; ++i)
1011 D(i) = abs(pow(dxs[i], 2) / det);
1017 : mesh(mesh_), length(length_)
1019 dim = mesh->Dimension();
1023 void CartesianPML::SetBoundaries()
1025 comp_dom_bdr.SetSize(dim, 2);
1026 dom_bdr.SetSize(dim, 2);
1028 mesh->GetBoundingBox(pmin, pmax);
1029 for (
int i = 0; i <
dim; i++)
1031 dom_bdr(i, 0) = pmin(i);
1032 dom_bdr(i, 1) = pmax(i);
1033 comp_dom_bdr(i, 0) = dom_bdr(i, 0) + length(i, 0);
1034 comp_dom_bdr(i, 1) = dom_bdr(i, 1) - length(i, 1);
1038 void CartesianPML::SetAttributes(
ParMesh *pmesh)
1041 for (
int i = 0; i < pmesh->
GetNBE(); ++i)
1046 int nrelem = pmesh->
GetNE();
1049 elems.SetSize(nrelem);
1052 for (
int i = 0; i < nrelem; ++i)
1055 bool in_pml =
false;
1062 int nrvert = vertices.
Size();
1065 for (
int iv = 0; iv < nrvert; ++iv)
1067 int vert_idx = vertices[iv];
1068 double *coords = pmesh->
GetVertex(vert_idx);
1069 for (
int comp = 0; comp <
dim; ++comp)
1071 if (coords[comp] > comp_dom_bdr(comp, 1) ||
1072 coords[comp] < comp_dom_bdr(comp, 0))
1088 void CartesianPML::StretchFunction(
const Vector &x,
1089 vector<complex<double>> &dxs)
1091 complex<double> zi = complex<double>(0., 1.);
1099 for (
int i = 0; i <
dim; ++i)
1104 coeff = n * c / k / pow(length(i, 1), n);
1105 dxs[i] = 1.0 + zi * coeff *
1110 coeff = n * c / k / pow(length(i, 0), n);
1111 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.
void SetOperator(const Operator &op)
Set the Operator and perform factorization.
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.
void SetPrintLevel(int print_lvl)
Set the error print level for MUMPS.
Pointer to an Operator of a specified type.
HYPRE_BigInt GlobalTrueVSize() const
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
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
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.
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 Mult(const Vector &x, Vector &y) const
Solve y = Op^{-1} x.
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).
void SetMatrixSymType(MatType mtype)
Set the matrix type.
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)
Jacobi smoothing for a given bilinear form (no matrix necessary).
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.
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)
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
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.
const Element * GetBdrElement(int i) const
void SetDiagonalBlock(int iblock, Operator *op)
Add a square block op in the block-entry (iblock, iblock).
MUMPS: A Parallel Sparse Direct Solver.
double f(const Vector &p)
bool Good() const
Return true if the command line options were parsed successfully.