39 return (0.25 * (2.0 + x[0]) - x[2]) * (x[2] + 0.25 * (2.0 + x[0]));
45 du[0] = 0.125 * (2.0 + x[0]) * x[1] * x[1];
46 du[1] = -0.125 * (2.0 + x[0]) * x[0] * x[1];
64 int main(
int argc,
char *argv[])
67 Mpi::Init(argc, argv);
68 int num_procs = Mpi::WorldSize();
69 int myid = Mpi::WorldRank();
76 int ser_ref_levels = 2;
77 int par_ref_levels = 1;
78 bool static_cond =
false;
79 bool visualization =
true;
82 args.
AddOption(&mesh_type,
"-mt",
"--mesh-type",
83 "Mesh type: 3 - Triangular, 4 - Quadrilateral.");
84 args.
AddOption(&mesh_order,
"-mo",
"--mesh-order",
85 "Geometric order of the curved mesh.");
86 args.
AddOption(&ser_ref_levels,
"-rs",
"--refine-serial",
87 "Number of times to refine the mesh uniformly in serial.");
88 args.
AddOption(&par_ref_levels,
"-rp",
"--refine-parallel",
89 "Number of times to refine the mesh uniformly in parallel.");
91 "Finite element order (polynomial degree) or -1 for"
92 " isoparametric space.");
93 args.
AddOption(&static_cond,
"-sc",
"--static-condensation",
"-no-sc",
94 "--no-static-condensation",
"Enable static condensation.");
95 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
97 "Enable or disable GLVis visualization.");
107 for (
int l = 0; l < ser_ref_levels; l++)
115 ParMesh pmesh(MPI_COMM_WORLD, *mesh);
117 for (
int l = 0; l < par_ref_levels; l++)
133 cout <<
"Number of unknowns: " << total_num_dofs << endl;
183 cout <<
"Size of linear system: "
206 if (myid == 0) { cout <<
"|u - u_h|_2 = " << error << endl; }
215 if (myid == 0) { cout <<
"|f - f_h|_2 = " << flux_err << endl; }
220 ostringstream mesh_name, sol_name, flux_name;
221 mesh_name <<
"mesh." << setfill(
'0') << setw(6) << myid;
222 sol_name <<
"sol." << setfill(
'0') << setw(6) << myid;
223 flux_name <<
"flux." << setfill(
'0') << setw(6) << myid;
225 ofstream mesh_ofs(mesh_name.str().c_str());
226 mesh_ofs.precision(8);
227 pmesh.
Print(mesh_ofs);
229 ofstream sol_ofs(sol_name.str().c_str());
230 sol_ofs.precision(8);
233 ofstream flux_ofs(flux_name.str().c_str());
234 flux_ofs.precision(8);
244 sol_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n";
245 sol_sock.precision(8);
246 sol_sock <<
"solution\n" << pmesh << x
247 <<
"window_title 'Solution'\n" << flush;
250 flux_sock <<
"parallel " << num_procs <<
" " << myid <<
"\n";
251 flux_sock.precision(8);
252 flux_sock <<
"solution\n" << pmesh << flux
254 <<
"window_geometry 402 0 400 350\n"
255 <<
"window_title 'Flux'\n" << flush;
269 mesh =
new Mesh(2, 12, 16, 8, 3);
312 mesh =
new Mesh(2, 8, 4, 8, 3);
339 MFEM_ABORT(
"Unrecognized mesh type " << type <<
"!");
354 if (fabs(x[1] + 1.0) < tol)
356 theta = 0.25 * M_PI * (x[0] - 2.0);
358 else if (fabs(x[0] - 1.0) < tol)
360 theta = 0.25 * M_PI * x[1];
362 else if (fabs(x[1] - 1.0) < tol)
364 theta = 0.25 * M_PI * (2.0 - x[0]);
366 else if (fabs(x[0] + 1.0) < tol)
368 theta = 0.25 * M_PI * (4.0 - x[1]);
372 cerr <<
"side not recognized "
373 << x[0] <<
" " << x[1] <<
" " << x[2] << endl;
378 r[2] = 0.25 * (2.0 * x[2] - 1.0) * (r[0] + 2.0);
385 double a = 17.0 - 2.0 * x[0] * (1.0 + x[0]);
386 s(0,0) = 0.5 + x[0] * x[0] * (8.0 / a - 0.5);
387 s(0,1) = x[0] * x[1] * (8.0 /
a - 0.5);
390 s(1,1) = 0.5 * x[0] * x[0] + 8.0 * x[1] * x[1] /
a;
int Size() const
Return the logical size of the array.
Class for domain integration L(v) := (f, v)
void duExact(const Vector &x, Vector &du)
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Conjugate gradient method.
void trans(const Vector &u, Vector &x)
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get()...
int AddQuad(int v1, int v2, int v3, int v4, int attr=1)
A coefficient that is constant across space and time.
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
void SetSize(int s)
Resize the vector to size s.
Pointer to an Operator of a specified type.
HYPRE_BigInt GlobalTrueVSize() const
void sigmaFunc(const Vector &x, DenseMatrix &s)
int AddTriangle(int v1, int v2, int v3, int attr=1)
Data type dense matrix using column-major storage.
void Transform(void(*f)(const Vector &, Vector &))
virtual void Save(std::ostream &out) const
Abstract parallel finite element space.
The BoomerAMG solver in hypre.
int AddVertex(double x, double y=0.0, double z=0.0)
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
A matrix coefficient with an optional scalar coefficient multiplier q. The matrix function can either...
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
void SetMaxIter(int max_it)
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
int AddBdrSegment(int v1, int v2, int attr=1)
FDualNumber< tbase > cos(const FDualNumber< tbase > &f)
cos([dual number])
A general vector function coefficient.
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
void SetRelTol(double rtol)
FDualNumber< tbase > sin(const FDualNumber< tbase > &f)
sin([dual number])
virtual void ComputeFlux(BilinearFormIntegrator &blfi, GridFunction &flux, bool wcoef=true, int subdomain=-1)
void FinalizeTopology(bool generate_bdr=true)
Finalize the construction of the secondary topology (connectivity) data of a Mesh.
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...
virtual double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL) const
double uExact(const Vector &x)
void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1) override
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
A general function coefficient.
void Mult(const double *x, double *y) const
Matrix vector multiplication.
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Arbitrary order H1-conforming (continuous) finite elements.
void Print(std::ostream &out=mfem::out) const override
Class for parallel grid function.
void fluxExact(const Vector &x, Vector &f)
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Wrapper for hypre's ParCSR matrix class.
Class for parallel meshes.
void ParseCheck(std::ostream &out=mfem::out)
double f(const Vector &p)
double sigma(const Vector &x)