33 return int(std::abs(1.0e5 * sin(seed * 1.1234 * M_PI)));
38int main(
int argc,
char *argv[])
48 const char *device_config =
"cpu";
49 bool visualization =
true;
52 bool deterministic =
true;
53 bool projectSolution =
false;
57 "Finite element order (polynomial degree) or -1 for"
58 " isoparametric space.");
59 args.
AddOption(&device_config,
"-d",
"--device",
60 "Device configuration string, see Device::Configure().");
61 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
63 "Enable or disable GLVis visualization.");
64 args.
AddOption(&numIter,
"-n",
"--num-iter",
"Number of hp-ref iterations");
65 args.
AddOption(&
dim,
"-dim",
"--dim",
"Mesh dimension (2 or 3)");
66 args.
AddOption(&deterministic,
"-det",
"--deterministic",
"-not-det",
67 "--not-deterministic",
68 "Whether to use deterministic random refinements");
69 args.
AddOption(&projectSolution,
"-proj",
"--project-solution",
"-no-proj",
71 "Whether to project a coefficient to solution");
89 Device device(device_config);
90 if (myid == 0) { device.
Print(); }
108 ParMesh pmesh(MPI_COMM_WORLD, mesh);
111 int par_ref_levels = 0;
112 for (
int l = 0; l < par_ref_levels; l++)
134 cout <<
"Using isoparametric FEs: " << fec->
Name() << endl;
143 const int fespaceDim = projectSolution ?
dim : 1;
152 const std::vector<char> hp_char = {
'h',
'p'};
154 for (
int iter=0; iter<numIter; ++iter)
156 const int r1 = deterministic ?
DetRand(seed) : rand();
157 const int r2 = deterministic ?
DetRand(seed) : rand();
158 const int elem = r1 % pmesh.
GetNE();
160 MPI_Bcast(&hp, 1, MPI_INT, 0, MPI_COMM_WORLD);
163 cout <<
"hp-refinement iteration " << iter <<
": "
164 << hp_char[hp] <<
"-refinement" << endl;
189 cout <<
"Number of finite element unknowns: " << size << endl;
190 cout <<
"Total number of h-refinements: " << numH
191 <<
"\nTotal number of p-refinements: " << numP
192 <<
"\nMaximum order " << maxP <<
"\n";
212 cout <<
"\n|| E_h - E ||_{L^2} = " << error <<
'\n' << endl;
256 a.FormLinearSystem(ess_tdof_list, x,
b, A, X, B);
279 cout << myid <<
": H1 continuity error " << h1error << endl;
280 MFEM_VERIFY(h1error < 1.0e-12,
"H1 continuity is not satisfied");
288 for (
int e=0; e<pmesh.
GetNE(); ++e)
293 MFEM_VERIFY(dofs.
Size() == 1,
"");
294 xo[dofs[0]] = p_elem;
301 ostringstream mesh_name, sol_name, order_name;
302 mesh_name <<
"mesh." << setfill(
'0') << setw(6) << myid;
303 sol_name <<
"sol." << setfill(
'0') << setw(6) << myid;
304 order_name <<
"order." << setfill(
'0') << setw(6) << myid;
306 ofstream mesh_ofs(mesh_name.str().c_str());
307 mesh_ofs.precision(8);
310 ofstream sol_ofs(sol_name.str().c_str());
311 sol_ofs.precision(8);
313 vis_x->Save(sol_ofs);
315 ofstream order_ofs(order_name.str().c_str());
316 order_ofs.precision(8);
326 sol_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n";
327 sol_sock.precision(8);
328 sol_sock <<
"solution\n" << pmesh << *vis_x << flush;
358 int Inf1, Inf2, NCFace;
365 auto &int_rule =
IntRules.
Get(FT->FaceGeom, 2 * faceOrder);
373 bool isNCSlave = FT->Elem2No >= 0 && NCFace >= 0;
374 bool isConforming = FT->Elem2No >= 0 && NCFace == -1;
375 if ((FT->Elem1No < FT->Elem2No && isConforming) || isNCSlave)
377 for (
int i = 0; i < nip; i++)
379 const auto &fip = int_rule.IntPoint(i);
382 FT->Loc1.Transform(fip, ip);
385 FT->Loc2.Transform(fip, ip);
388 const real_t err_i = std::abs(v1 - v2);
389 errorMax = std::max(errorMax, err_i);
400 if (!trueInterior) {
continue; }
404 const auto &int_rule =
IntRules.
Get(FT->FaceGeom, 2 * faceOrder);
407 for (
int i = 0; i < nip; i++)
409 const auto &fip = int_rule.IntPoint(i);
412 FT->Loc1.Transform(fip, ip);
415 FT->Loc2.Transform(fip, ip);
418 const real_t err_i = std::abs(v1 - v2);
419 errorMax = std::max(errorMax, err_i);
441 if (x.
Size() == 3) {
f(2) = 0.0; }
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
int Size() const
Return the logical size of the array.
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Conjugate gradient method.
void Mult(const Vector &b, Vector &x) const override
Iterative solution of the linear system using the Conjugate Gradient method.
void SetOperator(const Operator &op) override
Set/update the solver for the given operator.
A coefficient that is constant across space and time.
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Class for domain integration .
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
virtual const char * Name() const
int GetEdgeOrder(int edge, int variant=0) const
int GetFaceOrder(int face, int variant=0) const
Returns the polynomial degree of the i'th face finite element.
int GetElementOrder(int i) const
Returns the order of the i'th finite element.
Arbitrary order H1-conforming (continuous) finite elements.
The BoomerAMG solver in hypre.
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Class for integration point with weight.
int GetNPoints() const
Returns the number of the points in the integration rule.
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
void SetRelTol(real_t rtol)
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
void SetMaxIter(int max_it)
Arbitrary order "L2-conforming" discontinuous finite elements.
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
void GetFaceInfos(int Face, int *Inf1, int *Inf2) const
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
void Clear()
Clear the contents of the Mesh.
int GetNE() const
Returns number of elements.
int Dimension() const
Dimension of the reference space used within the elements.
bool FaceIsInterior(int FaceNo) const
Return true if the given face is interior.
static Mesh MakeCartesian3D(int nx, int ny, int nz, Element::Type type, real_t sx=1.0, real_t sy=1.0, real_t sz=1.0, bool sfc_ordering=true)
Creates a mesh for the parallelepiped [0,sx]x[0,sy]x[0,sz], divided into nx*ny*nz hexahedra if type =...
void GetNodes(Vector &node_coord) const
void EnsureNCMesh(bool simplices_nonconforming=false)
static Mesh MakeCartesian2D(int nx, int ny, Element::Type type, bool generate_edges=false, real_t sx=1.0, real_t sy=1.0, bool sfc_ordering=true)
Creates mesh for the rectangle [0,sx]x[0,sy], divided into nx*ny quadrilaterals if type = QUADRILATER...
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
static int WorldSize()
Return the size of MPI_COMM_WORLD.
static void Init(int &argc, char **&argv, int required=default_thread_required, int *provided=nullptr)
Singleton creation with Mpi::Init(argc, argv).
Pointer to an Operator of a specified type.
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
Reconstruct a solution vector x (e.g. a GridFunction) from the solution X of a constrained linear sys...
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
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...
bool Good() const
Return true if the command line options were parsed successfully.
Abstract parallel finite element space.
int GetMaxElementOrder() const override
Returns the maximum polynomial order over all elements globally.
void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1) const override
HYPRE_BigInt GlobalTrueVSize() const
int GetTrueVSize() const override
Return the number of local vector true dofs.
const Operator * GetProlongationMatrix() const override
void GetElementDofs(int i, Array< int > &dofs, DofTransformation &doftrans) const override
The same as GetElementDofs(), but with a user-allocated DofTransformation object. doftrans must be al...
void PRefineAndUpdate(const Array< pRefinement > &refs, bool want_transfer=true) override
const SparseMatrix * GetRestrictionMatrix() const override
Get the R matrix which restricts a local dof vector to true dof vector.
ParMesh * GetParMesh() const
void Update(bool want_transform=true) override
Class for parallel grid function.
void ExchangeFaceNbrData()
void Save(std::ostream &out) const override
real_t ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const override
Returns ||u_ex - u_h||_L2 in parallel for H1 or L2 elements.
real_t GetValue(int i, const IntegrationPoint &ip, int vdim=1) const override
ParFiniteElementSpace * ParFESpace() const
void ProjectCoefficient(Coefficient &coeff) override
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
std::unique_ptr< ParGridFunction > ProlongateToMaxOrder() const
Return a GridFunction with the values of this, prolongated to the maximum order of all elements in th...
Class for parallel meshes.
void ParPrint(std::ostream &out, const std::string &comments="") const
int GetNSharedFaces() const
Return the number of shared faces (3D), edges (2D), vertices (1D)
FaceElementTransformations * GetFaceElementTransformations(int FaceNo, int mask=31) override
bool FaceIsTrueInterior(int FaceNo) const
int GetSharedFace(int sface) const
Return the local face index for the given shared face.
FaceElementTransformations * GetSharedFaceTransformations(int sf, bool fill2=true)
Get the FaceElementTransformations for the given shared face (edge 2D) using the shared face index sf...
void Mult(const Vector &x, Vector &y) const override
Matrix vector multiplication.
A general vector function coefficient.
int Size() const
Returns the size of the vector.
void SetSize(int s)
Resize the vector to size s.
std::function< real_t(const Vector &)> f(real_t mass_coeff)
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
void f_exact(const Vector &x, Vector &f)
real_t CheckH1Continuity(ParGridFunction &x)