94int main(
int argc,
char *argv[])
97 const char *mesh_file =
"../../data/inline-quad.mesh";
101 bool visualization =
true;
104 bool static_cond =
false;
109 args.
AddOption(&mesh_file,
"-m",
"--mesh",
110 "Mesh file to use.");
112 "Finite element order (polynomial degree).");
113 args.
AddOption(&delta_order,
"-do",
"--delta-order",
114 "Order enrichment for DPG test space.");
116 "Epsilon coefficient");
117 args.
AddOption(&ref,
"-ref",
"--num-refinements",
118 "Number of uniform refinements");
119 args.
AddOption(&theta,
"-theta",
"--theta",
120 "Theta parameter for AMR");
121 args.
AddOption(&iprob,
"-prob",
"--problem",
"Problem case"
122 " 0: manufactured, 1: Erickson-Johnson ");
124 "Vector Coefficient beta");
125 args.
AddOption(&static_cond,
"-sc",
"--static-condensation",
"-no-sc",
126 "--no-static-condensation",
"Enable static condensation.");
127 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
128 "--no-visualization",
129 "Enable or disable GLVis visualization.");
130 args.
AddOption(&visport,
"-p",
"--send-port",
"Socket for GLVis.");
138 if (iprob > 1) { iprob = 1; }
141 if (
prob == prob_type::EJ)
143 mesh_file =
"../../data/inline-quad.mesh";
146 Mesh mesh(mesh_file, 1, 1);
148 MFEM_VERIFY(
dim > 1,
"Dimension = 1 is not supported in this example");
189 int test_order = order+delta_order;
213 trial_fes.
Append(sigma_fes);
214 trial_fes.
Append(hatu_fes);
215 trial_fes.
Append(hatf_fes);
226 a->StoreMatrices(
true);
230 TrialSpace::u_space, TestSpace::v_space);
234 TrialSpace::sigma_space, TestSpace::v_space);
238 TrialSpace::u_space, TestSpace::tau_space);
242 TrialSpace::sigma_space, TestSpace::tau_space);
246 TrialSpace::hatu_space, TestSpace::tau_space);
250 TrialSpace::hatf_space, TestSpace::v_space);
259 TestSpace::v_space, TestSpace::v_space);
262 TestSpace::v_space, TestSpace::v_space);
265 TestSpace::v_space, TestSpace::v_space);
268 TestSpace::tau_space, TestSpace::tau_space);
271 TestSpace::tau_space, TestSpace::tau_space);
290 std::cout <<
"\n Ref |"
295 <<
" Rate |" << std::endl;
296 std::cout << std::string(64,
'-') << std::endl;
298 if (static_cond) {
a->EnableStaticCondensation(); }
299 for (
int it = 0; it<=ref; it++)
311 ess_bdr_uhat = 1; ess_bdr_fhat = 0;
312 if (
prob == prob_type::EJ)
324 int n = ess_tdof_list_uhat.
Size();
325 int m = ess_tdof_list_fhat.
Size();
327 for (
int j = 0; j < n; j++)
329 ess_tdof_list[j] = ess_tdof_list_uhat[j]
333 for (
int j = 0; j < m; j++)
335 ess_tdof_list[j+n] = ess_tdof_list_fhat[j]
344 for (
int i = 0; i<trial_fes.
Size(); i++)
346 offsets[i+1] = trial_fes[i]->GetVSize();
347 dofs += trial_fes[i]->GetTrueVSize();
359 a->FormLinearSystem(ess_tdof_list,x,Ah,X,B);
378 a->RecoverFEMSolution(X,x);
386 real_t L2Error = sqrt(u_err*u_err + sigma_err*sigma_err);
388 Vector & residuals =
a->ComputeResidual(x);
391 real_t rate_err = (it) ?
dim*log(err0/L2Error)/log((
real_t)dof0/dofs) : 0.0;
392 real_t rate_res = (it) ?
dim*log(res0/residual)/log((
real_t)dof0/dofs) : 0.0;
398 std::ios oldState(
nullptr);
399 oldState.copyfmt(std::cout);
400 std::cout << std::right << std::setw(5) << it <<
" | "
401 << std::setw(10) << dof0 <<
" | "
402 << std::setprecision(3)
403 << std::setw(10) << std::scientific << err0 <<
" | "
404 << std::setprecision(2)
405 << std::setw(6) << std::fixed << rate_err <<
" | "
406 << std::setprecision(3)
407 << std::setw(10) << std::scientific << res0 <<
" | "
408 << std::setprecision(2)
409 << std::setw(6) << std::fixed << rate_res <<
" | "
411 std::cout.copyfmt(oldState);
415 const char * keys = (it == 0 &&
dim == 2) ?
"jRcm\n" :
nullptr;
418 "Numerical u", 0,0, 500, 500, keys);
420 "Numerical flux", 501,0,500, 500, keys);
430 for (
int iel = 0; iel<mesh.
GetNE(); iel++)
432 if (residuals[iel] > theta * max_resid)
434 elements_to_refine.
Append(iel);
439 for (
int i =0; i<trial_fes.
Size(); i++)
441 trial_fes[i]->Update(
false);
473 if (X.
Size() == 3) { z = X[2]; }
481 real_t denom = exp(-r2) - exp(-r1);
483 real_t g1 = exp(r2*(x-1.));
484 real_t g2 = exp(r1*(x-1.));
487 return g * cos(M_PI * y)/denom;
504 if (X.
Size() == 3) { z = X[2]; }
514 real_t denom = exp(-r2) - exp(-r1);
516 real_t g1 = exp(r2*(x-1.));
518 real_t g2 = exp(r1*(x-1.));
523 du[0] = g_x * cos(M_PI * y)/denom;
524 du[1] = -M_PI * g * sin(M_PI*y)/denom;
531 for (
int i = 0; i<du.
Size(); i++)
533 du[i] = M_PI * cos(
alpha);
545 if (X.
Size() == 3) { z = X[2]; }
554 real_t denom = exp(-r2) - exp(-r1);
556 real_t g1 = exp(r2*(x-1.));
559 real_t g2 = exp(r1*(x-1.));
563 real_t g_xx = g1_xx - g2_xx;
565 real_t u = g * cos(M_PI * y)/denom;
566 real_t u_xx = g_xx * cos(M_PI * y)/denom;
567 real_t u_yy = -M_PI * M_PI *
u;
575 return -M_PI*M_PI *
u * X.
Size();
599 for (
int i = 0; i<hatf.
Size(); i++)
613 for (
int i = 0; i<du.
Size(); i++)
615 s +=
beta[i] * du[i];
625 for (
int i = 0; i < mesh->
GetNE(); i++)
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
int Size() const
Return the logical size of the array.
void PartialSum()
Fill the entries of the array with the cumulative sum of the entries.
int Append(const T &el)
Append element 'el' to array, resize if necessary.
A class to handle Block diagonal preconditioners in a matrix-free implementation.
SparseMatrix & GetBlock(int i, int j)
Return a reference to block (i,j). Reference may be invalid if Aij(i,j) == NULL.
Array< int > & RowOffsets()
Return the row offsets for block starts.
int NumRowBlocks() const
Return the number of row blocks.
A class to handle Vectors in a block fashion.
Vector & GetBlock(int i)
Get the i-th vector in the block.
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.
Data type for scaled Jacobi-type smoother of sparse matrix.
Data type dense matrix using column-major storage.
for Raviart-Thomas elements
Class for domain integration .
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
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...
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 void Update(bool want_transform=true)
Reflect changes in the mesh: update number of DOFs, etc. Also, calculate GridFunction transformation ...
Mesh * GetMesh() const
Returns the mesh.
A general function coefficient.
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
Class for grid function - Vector with associated FE space.
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
virtual void MakeRef(FiniteElementSpace *f, real_t *v)
Make the GridFunction reference external data on a new FiniteElementSpace.
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()
void ProjectBdrCoefficientNormal(VectorCoefficient &vcoeff, const Array< int > &bdr_attr)
virtual void SetSpace(FiniteElementSpace *f)
Associate a new FiniteElementSpace with the GridFunction.
void ProjectBdrCoefficient(Coefficient &coeff, const Array< int > &attr)
Project a Coefficient on the GridFunction, modifying only DOFs on the boundary associated with the bo...
Arbitrary order H1-conforming (continuous) finite elements.
Arbitrary order "H^{1/2}-conforming" trace finite elements defined on the interface between mesh elem...
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.
A matrix coefficient that is constant in space and time.
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
int GetNE() const
Returns number of elements.
int Dimension() const
Dimension of the reference space used within the elements.
real_t GetElementVolume(int i)
Pointer to an Operator of a specified type.
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get().
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.
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Arbitrary order "H^{-1/2}-conforming" face finite elements defined on the interface between mesh elem...
Vector coefficient that is constant in space and time.
A general vector function coefficient.
void Neg()
(*this) = -(*this)
void SetSubVector(const Array< int > &dofs, const real_t value)
Set the entries listed in dofs to the given value.
real_t Norml2() const
Returns the l2 norm of the vector.
real_t Max() const
Returns the maximal element of the vector.
int Size() const
Returns the size of the vector.
void SetSize(int s)
Resize the vector to size s.
void setup_test_norm_coeffs(GridFunction &c1_gf, GridFunction &c2_gf)
real_t exact_hatu(const Vector &X)
real_t exact_laplacian_u(const Vector &X)
void exact_gradu(const Vector &X, Vector &du)
void exact_hatf(const Vector &X, Vector &hatf)
real_t exact_u(const Vector &X)
void exact_sigma(const Vector &X, Vector &sigma)
real_t f_exact(const Vector &X)
real_t sigma(const Vector &x)
void VisualizeField(socketstream &sock, const char *vishost, int visport, GridFunction &gf, const char *title, int x, int y, int w, int h, const char *keys, bool vec)
real_t u(const Vector &xvec)
void MultVVt(const Vector &v, DenseMatrix &vvt)
Make a matrix from a vector V.Vt.
std::function< real_t(const Vector &)> f(real_t mass_coeff)