69 #include "../common/mfem-common.hpp" 85 static const char *enum_str[] =
109 int main(
int argc,
char *argv[])
112 int myid = Mpi::WorldRank();
116 const char *mesh_file =
"../../data/inline-quad.mesh";
122 bool static_cond =
false;
125 bool visualization =
true;
126 bool paraview =
false;
129 args.
AddOption(&mesh_file,
"-m",
"--mesh",
130 "Mesh file to use.");
132 "Finite element order (polynomial degree).");
133 args.
AddOption(&delta_order,
"-do",
"--delta-order",
134 "Order enrichment for DPG test space.");
136 "Epsilon coefficient");
137 args.
AddOption(&ref,
"-ref",
"--num-refinements",
138 "Number of uniform refinements");
139 args.
AddOption(&theta,
"-theta",
"--theta",
140 "Theta parameter for AMR");
141 args.
AddOption(&iprob,
"-prob",
"--problem",
"Problem case" 142 " 0: lshape, 1: General");
144 "Vector Coefficient beta");
145 args.
AddOption(&static_cond,
"-sc",
"--static-condensation",
"-no-sc",
146 "--no-static-condensation",
"Enable static condensation.");
147 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
148 "--no-visualization",
149 "Enable or disable GLVis visualization.");
150 args.
AddOption(¶view,
"-paraview",
"--paraview",
"-no-paraview",
152 "Enable or disable ParaView visualization.");
163 if (iprob > 3) { iprob = 3; }
169 mesh_file =
"../../data/inline-quad.mesh";
172 Mesh mesh(mesh_file, 1, 1);
174 MFEM_VERIFY(
dim > 1,
"Dimension = 1 is not supported in this example");
210 ParMesh pmesh(MPI_COMM_WORLD, mesh);
244 int test_order = order+delta_order;
266 trial_fes.
Append(sigma_fes);
267 trial_fes.
Append(hatu_fes);
268 trial_fes.
Append(hatf_fes);
273 a->StoreMatrices(
true);
277 TrialSpace::u_space, TestSpace::v_space);
281 TrialSpace::sigma_space, TestSpace::v_space);
285 TrialSpace::u_space, TestSpace::tau_space);
289 TrialSpace::sigma_space, TestSpace::tau_space);
293 TrialSpace::hatu_space, TestSpace::tau_space);
297 TrialSpace::hatf_space, TestSpace::v_space);
312 TestSpace::v_space, TestSpace::v_space);
315 TestSpace::v_space, TestSpace::v_space);
318 TestSpace::v_space, TestSpace::v_space);
321 TestSpace::tau_space, TestSpace::tau_space);
324 TestSpace::tau_space, TestSpace::tau_space);
350 std::cout <<
" Ref |" 357 std::cout <<
" Residual |" 359 <<
" CG it |" << endl;
360 std::cout << std::string((
exact_known) ? 72 : 50,
'-')
364 if (static_cond) {
a->EnableStaticCondensation(); }
384 for (
int it = 0; it<=ref; it++)
415 int n = ess_tdof_list_uhat.
Size();
416 int m = ess_tdof_list_fhat.
Size();
418 for (
int j = 0; j < n; j++)
420 ess_tdof_list[j] = ess_tdof_list_uhat[j]
424 for (
int j = 0; j < m; j++)
426 ess_tdof_list[j+n] = ess_tdof_list_fhat[j]
450 a->FormLinearSystem(ess_tdof_list,x,Ah,X,B);
462 amg1->SetPrintLevel(0);
463 M.SetDiagonalBlock(0,amg0);
464 M.SetDiagonalBlock(1,amg1);
470 M.SetDiagonalBlock(skip,amg2);
483 M.SetDiagonalBlock(skip+1,prec);
494 a->RecoverFEMSolution(X,x);
495 Vector & residuals =
a->ComputeResidual(x);
497 double residual = residuals.
Norml2();
498 double maxresidual = residuals.
Max();
500 double gresidual = residual * residual;
502 MPI_Allreduce(MPI_IN_PLACE,&maxresidual,1,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD);
503 MPI_Allreduce(MPI_IN_PLACE,&gresidual,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
505 gresidual = sqrt(gresidual);
508 for (
int iel = 0; iel<pmesh.
GetNE(); iel++)
510 if (residuals[iel] > theta * maxresidual)
512 elements_to_refine.
Append(iel);
524 double L2Error = 0.0;
525 double rate_err = 0.0;
530 L2Error = sqrt(u_err*u_err + sigma_err*sigma_err);
531 rate_err = (it) ?
dim*log(err0/L2Error)/log((
double)dof0/dofs) : 0.0;
534 double rate_res = (it) ?
dim*log(res0/gresidual)/log((
double)dof0/dofs) : 0.0;
541 std::ios oldState(
nullptr);
542 oldState.copyfmt(std::cout);
543 std::cout << std::right << std::setw(5) << it <<
" | " 544 << std::setw(10) << dof0 <<
" | ";
547 std::cout << std::setprecision(3) << std::setw(10)
548 << std::scientific << err0 <<
" | " 549 << std::setprecision(2)
550 << std::setw(6) << std::fixed << rate_err <<
" | " ;
552 std::cout << std::setprecision(3)
553 << std::setw(10) << std::scientific << res0 <<
" | " 554 << std::setprecision(2)
555 << std::setw(6) << std::fixed << rate_res <<
" | " 556 << std::setw(6) << std::fixed << num_iter <<
" | " 558 std::cout.copyfmt(oldState);
563 const char * keys = (it == 0 &&
dim == 2) ?
"cgRjmlk\n" :
nullptr;
567 "Numerical u", 0,0, 500, 500, keys);
569 "Numerical flux", 501,0,500, 500, keys);
575 paraview_dc->
SetTime((
double)it);
585 for (
int i =0; i<trial_fes.
Size(); i++)
587 trial_fes[i]->Update(
false);
624 if (X.
Size() == 3) { z = X[2]; }
629 double alpha = M_PI * (x + y + z);
638 double denom = exp(-r2) - exp(-r1);
640 double g1 = exp(r2*(x-1.));
641 double g2 = exp(r1*(x-1.));
643 return g * cos(M_PI * y)/denom;
648 double r = sqrt(x*x+y*y);
653 MFEM_ABORT(
"Wrong code path");
664 if (X.
Size() == 3) { z = X[2]; }
671 double alpha = M_PI * (x + y + z);
672 for (
int i = 0; i<du.
Size(); i++)
674 du[i] = M_PI * cos(
alpha);
683 double denom = exp(-r2) - exp(-r1);
685 double g1 = exp(r2*(x-1.));
687 double g2 = exp(r1*(x-1.));
690 double g_x = g1_x - g2_x;
692 double u_x = g_x * cos(M_PI * y)/denom;
693 double u_y = -M_PI * g * sin(M_PI*y)/denom;
700 double r = sqrt(x*x+y*y);
702 double denom = r*
alpha;
708 MFEM_ABORT(
"Wrong code path");
718 if (X.
Size() == 3) { z = X[2]; }
723 double alpha = M_PI * (x + y + z);
725 return - M_PI*M_PI *
u * X.
Size();
733 double denom = exp(-r2) - exp(-r1);
735 double g1 = exp(r2*(x-1.));
737 double g1_xx = r2*g1_x;
738 double g2 = exp(r1*(x-1.));
740 double g2_xx = r1*g2_x;
742 double g_xx = g1_xx - g2_xx;
744 double u = g * cos(M_PI * y)/denom;
745 double u_xx = g_xx * cos(M_PI * y)/denom;
746 double u_yy = -M_PI * M_PI *
u;
752 double r = sqrt(x*x+y*y);
758 MFEM_ABORT(
"Wrong code path");
784 for (
int i = 0; i<hatf.
Size(); i++)
786 hatf[i] = beta_val[i] *
u -
sigma[i];
801 for (
int i = 0; i<du.
Size(); i++)
803 s += beta_val[i] * du[i];
841 beta_val(0) = exp(x)*sin(y);
842 beta_val(1) = exp(x)*cos(y);
855 for (
int i = 0; i < pmesh->
GetNE(); i++)
858 double c1 = min(
epsilon/volume, 1.);
859 double c2 = min(1./
epsilon, 1./volume);
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Class for domain integration L(v) := (f, v)
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Conjugate gradient method.
The Auxiliary-space Maxwell Solver in hypre.
double exact_u(const Vector &X)
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
The Auxiliary-space Divergence Solver in hypre.
void exact_sigma(const Vector &X, Vector &sigma)
A class to handle Vectors in a block fashion.
void SetDataFormat(VTKFormat fmt)
A coefficient that is constant across space and time.
void beta_function(const Vector &X, Vector &beta_val)
Array< int > & RowOffsets()
Return the row offsets for block starts.
void PrintOptions(std::ostream &out) const
Print the options.
void setup_test_norm_coeffs(ParGridFunction &c1_gf, ParGridFunction &c2_gf)
int Dimension() const
Dimension of the reference space used within the elements.
void SetSize(int s)
Resize the vector to size s.
Helper class for ParaView visualization data.
void PrintUsage(std::ostream &out) const
Print the usage message.
virtual void Update(bool want_transform=true)
Pointer to an Operator of a specified type.
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
int Size() const
Returns the size of the vector.
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
virtual void SetSpace(FiniteElementSpace *f)
Associate a new FiniteElementSpace with the ParGridFunction.
bool Good() const
Return true if the command line options were parsed successfully.
Abstract parallel finite element space.
(Q div u, div v) for RT elements
void exact_gradu(const Vector &X, Vector &du)
Operator & GetBlock(int i, int j)
Return a reference to block i,j.
double bdr_data(const Vector &X)
Matrix coefficient defined as the outer product of two vector coefficients.
virtual void MakeRef(FiniteElementSpace *f, double *v)
Make the ParGridFunction reference external data on a new FiniteElementSpace.
The BoomerAMG solver in hypre.
void ProjectBdrCoefficientNormal(VectorCoefficient &vcoeff, Array< int > &bdr_attr)
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
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. ...
int Append(const T &el)
Append element 'el' to array, resize if necessary.
void EnsureNCMesh(bool simplices_nonconforming=false)
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
int GetNumIterations() const
Returns the number of iterations taken during the last call to Mult()
Arbitrary order "H^{1/2}-conforming" trace finite elements defined on the interface between mesh elem...
void SetPrintLevel(int print_level)
ParMesh * GetParMesh() const
void SetMaxIter(int max_it)
void SetHighOrderOutput(bool high_order_output_)
void VisualizeField(socketstream &sock, const char *vishost, int visport, ParGridFunction &gf, const char *title, int x=0, int y=0, int w=400, int h=400, bool vec=false)
void SetTime(double t)
Set physical time (for time-dependent simulations)
virtual double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const
HYPRE_BigInt GlobalTrueVSize() const
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
virtual DofTransformation * GetElementDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom in array dofs for i'th element.
A general vector function coefficient.
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
void SetRelTol(double rtol)
void ProjectBdrCoefficient(Coefficient *coeff[], VectorCoefficient *vcoeff, Array< int > &attr)
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Vector coefficient defined as a product of scalar and vector coefficients.
ParFiniteElementSpace * ParFESpace() const
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.
double Max() const
Returns the maximal element of the vector.
int GetNE() const
Returns number of elements.
virtual int GetTrueVSize() const
Return the number of local vector true dofs.
int main(int argc, char *argv[])
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get()...
double exact_laplacian_u(const Vector &X)
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
double Norml2() const
Returns the l2 norm of the vector.
Abstract class for hypre's solvers and preconditioners.
int Size() const
Return the logical size of the array.
void SetLevelsOfDetail(int levels_of_detail_)
void Clear()
Clear the contents of the Mesh.
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
A general function coefficient.
Arbitrary order "H^{-1/2}-conforming" face finite elements defined on the interface between mesh elem...
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Arbitrary order H1-conforming (continuous) finite elements.
double sigma(const Vector &x)
virtual void Save() override
double u(const Vector &xvec)
Class for parallel grid function.
double exact_hatu(const Vector &X)
Wrapper for hypre's ParCSR matrix class.
A class to handle Block systems in a matrix-free implementation.
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Class for parallel meshes.
double GetElementVolume(int i)
void SetPrefixPath(const std::string &prefix)
Set the path where the DataCollection will be saved.
void exact_hatf(const Vector &X, Vector &hatf)
Vector & GetBlock(int i)
Get the i-th vector in the block.
Arbitrary order "L2-conforming" discontinuous finite elements.
double f(const Vector &p)
double f_exact(const Vector &X)