33 return int(std::abs(1.0e5 * sin(seed * 1.1234 * M_PI)));
38int main(
int argc,
char *argv[])
42 const char *device_config =
"cpu";
43 bool visualization =
true;
46 bool deterministic =
true;
47 bool projectSolution =
false;
51 "Finite element order (polynomial degree) or -1 for"
52 " isoparametric space.");
53 args.
AddOption(&device_config,
"-d",
"--device",
54 "Device configuration string, see Device::Configure().");
55 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
57 "Enable or disable GLVis visualization.");
58 args.
AddOption(&numIter,
"-n",
"--num-iter",
"Number of hp-ref iterations");
59 args.
AddOption(&
dim,
"-dim",
"--dim",
"Mesh dimension (2 or 3)");
60 args.
AddOption(&deterministic,
"-det",
"--deterministic",
"-not-det",
61 "--not-deterministic",
62 "Whether to use deterministic random refinements");
63 args.
AddOption(&projectSolution,
"-proj",
"--project-solution",
"-no-proj",
65 "Whether to project a coefficient to solution");
76 Device device(device_config);
106 cout <<
"Using isoparametric FEs: " << fec->
Name() << endl;
114 const int fespaceDim = projectSolution ?
dim : 1;
122 const std::vector<char> hp_char = {
'h',
'p'};
124 for (
int iter=0; iter<numIter; ++iter)
126 const int r1 = deterministic ?
DetRand(seed) : rand();
127 const int r2 = deterministic ?
DetRand(seed) : rand();
128 const int elem = r1 % mesh.
GetNE();
129 const int hp = r2 % 2;
131 cout <<
"hp-refinement iteration " << iter <<
": "
132 << hp_char[hp] <<
"-refinement" << endl;
154 cout <<
"Number of finite element unknowns: " << size << endl;
157 cout <<
"Total number of h-refinements: " << numH
158 <<
"\nTotal number of p-refinements: " << numP
159 <<
"\nMaximum order " << maxP <<
"\n";
176 cout <<
"\n|| E_h - E ||_{L^2} = " << error <<
'\n' << endl;
219 a.FormLinearSystem(ess_tdof_list, x,
b, A, X, B);
223#ifndef MFEM_USE_SUITESPARSE
226 PCG(*A, M, B, X, 1, 200, 1e-12, 0.0);
230 umf_solver.
Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
232 umf_solver.
Mult(B, X);
237 a.RecoverFEMSolution(X,
b, x);
243 cout <<
"H1 continuity error " << h1error << endl;
244 MFEM_VERIFY(h1error < 1.0e-12,
"H1 continuity is not satisfied");
252 for (
int e=0; e<mesh.
GetNE(); ++e)
257 MFEM_VERIFY(dofs.
Size() == 1,
"");
258 xo[dofs[0]] = p_elem;
264 ofstream mesh_ofs(
"refined.mesh");
265 mesh_ofs.precision(8);
266 mesh.
Print(mesh_ofs);
267 ofstream sol_ofs(
"sol.gf");
268 sol_ofs.precision(8);
269 vis_x->Save(sol_ofs);
270 ofstream order_ofs(
"order.gf");
271 order_ofs.precision(8);
280 sol_sock.precision(8);
281 sol_sock <<
"solution\n" << mesh << *vis_x << flush;
307 int Inf1, Inf2, NCFace;
314 auto &int_rule =
IntRules.
Get(FT->FaceGeom, 2 * faceOrder);
322 bool isNCSlave = FT->Elem2No >= 0 && NCFace >= 0;
323 bool isConforming = FT->Elem2No >= 0 && NCFace == -1;
324 if ((FT->Elem1No < FT->Elem2No && isConforming) || isNCSlave)
326 for (
int i = 0; i < nip; i++)
328 const auto &fip = int_rule.IntPoint(i);
331 FT->Loc1.Transform(fip, ip);
334 FT->Loc2.Transform(fip, ip);
337 const real_t err_i = std::abs(v1 - v2);
338 errorMax = std::max(errorMax, err_i);
362 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.
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
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
DofTransformation * GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element 'elem'. The returned indices are offsets into an ldo...
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.
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1) const
Get a list of essential true dofs, ess_tdof_list, corresponding to the boundary attributes marked in ...
virtual const Operator * GetProlongationMatrix() const
virtual void Update(bool want_transform=true)
Reflect changes in the mesh: update number of DOFs, etc. Also, calculate GridFunction transformation ...
int GetElementOrder(int i) const
Returns the order of the i'th finite element.
Mesh * GetMesh() const
Returns the mesh.
virtual void PRefineAndUpdate(const Array< pRefinement > &refs, bool want_transfer=true)
virtual const SparseMatrix * GetHpRestrictionMatrix() const
The returned SparseMatrix is owned by the FiniteElementSpace.
virtual int GetMaxElementOrder() const
Return the maximum polynomial order over all elements.
Data type for Gauss-Seidel smoother of sparse matrix.
Class for grid function - Vector with associated FE space.
virtual real_t GetValue(int i, const IntegrationPoint &ip, int vdim=1) const
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
virtual real_t ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const
Returns ||exsol - u_h||_L2 for scalar or vector H1 or L2 elements.
FiniteElementSpace * FESpace()
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
std::unique_ptr< GridFunction > ProlongateToMaxOrder() const
Return a GridFunction with the values of this, prolongated to the maximum order of all elements in th...
Arbitrary order H1-conforming (continuous) finite elements.
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.
Arbitrary order "L2-conforming" discontinuous finite elements.
virtual FaceElementTransformations * GetFaceElementTransformations(int FaceNo, int mask=31)
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)
virtual void Print(std::ostream &os=mfem::out, const std::string &comments="") const
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...
Pointer to an Operator of a specified type.
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
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.
void Mult(const Vector &x, Vector &y) const override
Matrix vector multiplication.
Direct sparse solver using UMFPACK.
real_t Control[UMFPACK_CONTROL]
void SetOperator(const Operator &op) override
Factorize the given Operator op which must be a SparseMatrix.
void Mult(const Vector &b, Vector &x) const override
Direct solution of the linear system using UMFPACK.
A general vector function coefficient.
int Size() const
Returns the size of the vector.
void SetSize(int s)
Resize the vector to size s.
real_t CheckH1Continuity(GridFunction &x)
void f_exact(const Vector &x, Vector &f)
void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x, int print_iter, int max_num_iter, real_t RTOLERANCE, real_t ATOLERANCE)
Preconditioned conjugate gradient method. (tolerances are squared)
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)