67 int main(
int argc,
char *argv[])
72 const char *mesh_file =
"../data/star-hilbert.mesh";
75 double max_elem_error = 5.0e-3;
76 double hysteresis = 0.15;
79 bool visualization =
true;
83 args.
AddOption(&mesh_file,
"-m",
"--mesh",
86 "Problem setup to use: 0 = spherical front, 1 = ball.");
88 "Number of solution features (fronts/balls).");
90 "Finite element order (polynomial degree).");
91 args.
AddOption(&max_elem_error,
"-e",
"--max-err",
92 "Maximum element error");
93 args.
AddOption(&hysteresis,
"-y",
"--hysteresis",
94 "Derefinement safety coefficient.");
95 args.
AddOption(&ref_levels,
"-r",
"--ref-levels",
96 "Number of initial uniform refinement levels.");
97 args.
AddOption(&nc_limit,
"-l",
"--nc-limit",
98 "Maximum level of hanging nodes.");
99 args.
AddOption(&t_final,
"-tf",
"--t-final",
100 "Final time; start time is 0.");
101 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
102 "--no-visualization",
103 "Enable or disable GLVis visualization.");
104 args.
AddOption(&visit,
"-visit",
"--visit-datafiles",
"-no-visit",
105 "--no-visit-datafiles",
106 "Save data files for VisIt (visit.llnl.gov) visualization.");
118 Mesh mesh(mesh_file, 1, 1);
127 if (ref_levels > 0) { ref_levels--; }
131 for (
int l = 0; l < ref_levels; l++)
140 "Boundary attributes required in the mesh.");
168 char vishost[] =
"localhost";
174 sout.
open(vishost, visport);
177 cout <<
"Unable to connect to GLVis server at "
178 << vishost <<
':' << visport << endl;
179 cout <<
"GLVis visualization disabled.\n";
180 visualization =
false;
219 for (
double time = 0.0; time < t_final + 1e-10; time += 0.01)
221 cout <<
"\nTime " << time <<
"\n\nRefinement:" << endl;
233 for (
int ref_it = 1; ; ref_it++)
235 cout <<
"Iteration: " << ref_it <<
", number of unknowns: "
254 #ifndef MFEM_USE_SUITESPARSE
256 PCG(A, M, B, X, 0, 500, 1e-12, 0.0);
259 umf_solver.
Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
261 umf_solver.
Mult(B, X);
272 sout <<
"solution\n" << mesh << x << flush;
299 if (derefiner.
Apply(mesh))
301 cout <<
"\nDerefined elements." << endl;
339 double front(
double x,
double y,
double z,
double t,
int)
341 double r = sqrt(x*x + y*y + z*z);
342 return exp(-0.5*pow((r - t)/
alpha, 2));
347 double x2 = x*x, y2 = y*y, z2 = z*z, t2 = t*t;
348 double r = sqrt(x2 + y2 + z2);
350 return -exp(-0.5*pow((r - t)/alpha, 2)) / a4 *
351 (-2*t*(x2 + y2 + z2 - (dim-1)*a2/2)/r + x2 + y2 + z2 + t2 - dim*a2);
355 double ball(
double x,
double y,
double z,
double t,
int)
357 double r = sqrt(x*x + y*y + z*z);
358 return -atan(2*(r - t)/
alpha);
363 double x2 = x*x, y2 = y*y, z2 = z*z, t2 = 4*t*t;
364 double r = sqrt(x2 + y2 + z2);
366 double den = pow(-a2 - 4*(x2 + y2 + z2 - 2*r*t) - t2, 2.0);
367 return (dim == 2) ? 2*alpha*(a2 + t2 - 4*x2 - 4*y2)/r/den
368 : 4*alpha*(a2 + t2 - 4*r*t)/r/den;
372 template<
typename F0,
typename F1>
376 double x = pt(0), y = pt(1), z = 0.0;
377 if (dim == 3) { z = pt(2); }
383 return f0(x, y, z, t, dim);
390 double x0 = 0.5*cos(2*M_PI * i / nfeatures);
391 double y0 = 0.5*sin(2*M_PI * i / nfeatures);
392 sum += f0(x - x0, y - y0, z, t, dim);
402 double x0 = 0.5*cos(2*M_PI * i / nfeatures + M_PI*t);
403 double y0 = 0.5*sin(2*M_PI * i / nfeatures + M_PI*t);
404 sum += f1(x - x0, y - y0, z, 0.25, dim);
double front_laplace(double x, double y, double z, double t, int dim)
int Size() const
Logical size of the array.
Class for domain integration L(v) := (f, v)
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
double bdr_func(const Vector &pt, double t)
Class for grid function - Vector with associated FE space.
virtual void Update(bool want_transform=true)
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
Subclass constant coefficient.
void SetNCLimit(int nc_limit)
Set the maximum ratio of refinement levels of adjacent elements (0 = unlimited).
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)
virtual void Save()
Save the collection and a VisIt root file.
bool Stop() const
Check if STOP action is requested, e.g. stopping criterion is satisfied.
int main(int argc, char *argv[])
double rhs_func(const Vector &pt, double t)
void UpdateProblem(Mesh &mesh, FiniteElementSpace &fespace, GridFunction &x, BilinearForm &a, LinearForm &b)
Data type for Gauss-Seidel smoother of sparse matrix.
Direct sparse solver using UMFPACK.
bool Apply(Mesh &mesh)
Perform the mesh operation.
virtual void Reset()
Reset the associated estimator.
void EnsureNCMesh(bool simplices_nonconforming=false)
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
virtual void Reset()
Reset the associated estimator.
Mesh refinement operator using an error threshold.
Data collection with VisIt I/O routines.
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
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)
void PrintUsage(std::ostream &out) const
void SetTime(double t)
Set physical time (for time-dependent simulations)
void SetNCLimit(int nc_limit)
Set the maximum ratio of refinement levels of adjacent elements (0 = unlimited).
int SpaceDimension() const
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
double ball_laplace(double x, double y, double z, double t, int dim)
double Control[UMFPACK_CONTROL]
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
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)
void PreferConformingRefinement()
Use conforming refinement, if possible (triangles, tetrahedra) – this is the default.
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection and update the root file.
NURBSExtension * NURBSext
Optional NURBS mesh extension.
virtual void Finalize(bool refine=false, bool fix_orientation=false)
Finalize the construction of a general Mesh.
void PrintOptions(std::ostream &out) const
double ball(double x, double y, double z, double t, int)
int open(const char hostname[], int port)
void SetThreshold(double thresh)
Set the de-refinement threshold. The default value is zero.
The ZienkiewiczZhuEstimator class implements the Zienkiewicz-Zhu error estimation procedure...
void SetLocalErrorGoal(double err_goal)
Set the local stopping criterion: stop when local_err_i <= local_err_goal. The default value is zero...
class for C-function coefficient
virtual void UpdatesFinished()
Free the GridFunction update operator (if any), to save memory.
Arbitrary order H1-conforming (continuous) finite elements.
void SetTotalErrorFraction(double fraction)
Set the total fraction used in the computation of the threshold. The default value is 1/2...
double front(double x, double y, double z, double t, int)
double composite_func(const Vector &pt, double t, F0 f0, F1 f1)
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
void ProjectBdrCoefficient(Coefficient &coeff, Array< int > &attr)
Project a Coefficient on the GridFunction, modifying only DOFs on the boundary associated with the bo...
De-refinement operator using an error threshold.
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a SparseMatrix.