87 int main(
int argc,
char *argv[])
90 const char *mesh_file =
"../data/star.mesh";
94 bool visualization =
true;
95 bool verification =
false;
98 args.
AddOption(&mesh_file,
"-m",
"--mesh",
101 "Finite element order (polynomial degree) or -1 for"
102 " isoparametric space.");
103 args.
AddOption(&num_refs,
"-r",
"--refs",
104 "Number of uniform refinements");
105 args.
AddOption(&alpha,
"-alpha",
"--alpha",
106 "Fractional exponent");
107 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
108 "--no-visualization",
109 "Enable or disable GLVis visualization.");
110 args.
AddOption(&verification,
"-ver",
"--verification",
"-no-ver",
112 "Use sinusoidal function (f) for analytic comparison.");
122 int progress_steps = 1;
126 const int power_of_laplace = floor(alpha);
127 double exponent_to_approximate = alpha - power_of_laplace;
128 bool integer_order =
false;
130 if (abs(exponent_to_approximate) > 1e-12)
132 mfem::out <<
"Approximating the fractional exponent "
133 << exponent_to_approximate
140 alpha = exponent_to_approximate + power_of_laplace;
144 integer_order =
true;
145 mfem::out <<
"Treating integer order PDE." << endl;
149 Mesh mesh(mesh_file, 1, 1);
153 for (
int i = 0; i < num_refs; i++)
161 cout <<
"Number of finite element unknowns: "
177 for (
int i=0; i<x.Size(); i++)
179 val *= sin(M_PI*x(i));
181 return pow(x.Size()*pow(M_PI,2),
alpha) * val;
217 if (power_of_laplace > 0)
238 mfem::out <<
"\nComputing (-Δ) ^ -" << power_of_laplace
239 <<
" ( f ) " << endl;
240 for (
int i = 0; i < power_of_laplace; i++)
243 PCG(*Op, M, B, X, 3, 300, 1e-12, 0.0);
246 if (i == power_of_laplace - 1)
250 if (integer_order && verification)
259 fout.
open(vishost, visport);
261 oss_f.str(
""); oss_f.clear();
262 oss_f <<
"Step " << progress_steps++ <<
": Solution of PDE -Δ ^ "
265 fout <<
"solution\n" << mesh << g
266 <<
"window_title '" << oss_f.str() <<
"'" << flush;
296 ostringstream oss_x, oss_u;
299 xout.
open(vishost, visport);
301 uout.
open(vishost, visport);
306 for (
int i = 0; i < coeffs.
Size(); i++)
308 mfem::out <<
"\nSolving PDE -Δ u + " << -poles[i]
309 <<
" u = " << coeffs[i] <<
" g " << endl;
330 PCG(*A, M, B, X, 3, 300, 1e-12, 0.0);
342 oss_x.str(
""); oss_x.clear();
343 oss_x <<
"Step " << progress_steps
344 <<
": Solution of PDE -Δ u + " << -poles[i]
345 <<
" u = " << coeffs[i] <<
" g";
346 xout <<
"solution\n" << mesh << x
347 <<
"window_title '" << oss_x.str() <<
"'" << flush;
349 oss_u.str(
""); oss_u.clear();
350 oss_u <<
"Step " << progress_steps + 1
351 <<
": Solution of fractional PDE (-Δ)^" << alpha
353 uout <<
"solution\n" << mesh << u
354 <<
"window_title '" << oss_u.str() <<
"'"
365 auto solution = [] (
const Vector &x)
368 for (
int i=0; i<x.
Size(); i++)
370 val *= sin(M_PI*x(i));
381 analytic_solution =
"sin(π x)";
382 expected_mesh =
"inline_segment.mesh";
385 analytic_solution =
"sin(π x) sin(π y)";
386 expected_mesh =
"inline_quad.mesh";
389 analytic_solution =
"sin(π x) sin(π y) sin(π z)";
390 expected_mesh =
"inline_hex.mesh";
395 <<
"\n\nSolution Verification in "<< dim <<
"D \n\n"
396 <<
"Analytic solution : " << analytic_solution <<
"\n"
397 <<
"Expected mesh : " << expected_mesh <<
"\n"
398 <<
"Your mesh : " << mesh_file <<
"\n"
399 <<
"L2 error : " << l2_error <<
"\n\n"
400 << string(80,
'=') << endl;
int Size() const
Return the logical size of the array.
Class for domain integration L(v) := (f, v)
Class for grid function - Vector with associated FE space.
A coefficient that is constant across space and time.
Pointer to an Operator of a specified type.
int Size() const
Returns the size of the vector.
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Get a list of essential true dofs, ess_tdof_list, corresponding to the boundary attributes marked in ...
Data type for Gauss-Seidel smoother of sparse matrix.
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
double analytic_solution(const Vector &x)
void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x, int print_iter, int max_num_iter, double RTOLERANCE, double ATOLERANCE)
Preconditioned conjugate gradient method. (tolerances are squared)
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
void PrintUsage(std::ostream &out) const
Print the usage message.
void SetSubVectorComplement(const Array< int > &dofs, const double val)
Set all vector entries NOT in the dofs Array to the given val.
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
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 ComputePartialFractionApproximation(double &alpha, Array< double > &coeffs, Array< double > &poles, double lmax=1000., double tol=1e-10, int npoints=1000, int max_order=100)
virtual double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const
void MultTranspose(const Vector &x, Vector &y) const
Multiply a vector with the transposed matrix. y = At * x.
void PrintOptions(std::ostream &out) const
Print the options.
int open(const char hostname[], int port)
Open the socket stream on 'port' at 'hostname'.
A general function coefficient.
Arbitrary order H1-conforming (continuous) finite elements.
double u(const Vector &xvec)
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
virtual const SparseMatrix * GetRestrictionMatrix() const
The returned SparseMatrix is owned by the FiniteElementSpace.
double f(const Vector &p)
bool Good() const
Return true if the command line options were parsed successfully.