45 int main(
int argc,
char *argv[])
48 const char *mesh_file =
"../data/inline-quad.mesh";
51 bool visualization = 1;
54 args.
AddOption(&mesh_file,
"-m",
"--mesh",
56 args.
AddOption(&ref_levels,
"-r",
"--refine",
57 "Number of times to refine the mesh uniformly.");
59 "Finite element order (polynomial degree).");
60 args.
AddOption(&
freq,
"-f",
"--frequency",
"Set the frequency for the exact"
62 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
64 "Enable or disable GLVis visualization.");
71 Mesh mesh(mesh_file, 1, 1);
77 for (
int lev = 0; lev < ref_levels; lev++)
100 cout <<
"Number of H(Curl) unknowns: " << size << endl;
137 sigmaMat(0,0) = 2.0; sigmaMat(1,1) = 2.0; sigmaMat(2,2) = 2.0;
138 sigmaMat(0,2) = 0.0; sigmaMat(2,0) = 0.0;
139 sigmaMat(0,1) = M_SQRT1_2; sigmaMat(1,0) = M_SQRT1_2;
140 sigmaMat(1,2) = M_SQRT1_2; sigmaMat(2,1) = M_SQRT1_2;
161 #ifndef MFEM_USE_SUITESPARSE
165 PCG(*A, M, B, X, 1, 500, 1e-12, 0.0);
170 umf_solver.
Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
172 umf_solver.
Mult(B, X);
181 cout <<
"\n|| E_h - E ||_{H(Curl)} = " << error <<
'\n' << endl;
188 ofstream mesh_ofs(
"refined.mesh");
189 mesh_ofs.precision(8);
190 mesh.
Print(mesh_ofs);
191 ofstream sol_ofs(
"sol.gf");
192 sol_ofs.precision(8);
215 dy_sock.precision(8);
216 dz_sock.precision(8);
218 Vector xVec(3); xVec = 0.0; xVec(0) = 1;
219 Vector yVec(3); yVec = 0.0; yVec(1) = 1;
220 Vector zVec(3); zVec = 0.0; zVec(2) = 1;
246 x_sock <<
"solution\n" << mesh << xComp << flush
247 <<
"window_title 'X component'" << endl;
248 y_sock <<
"solution\n" << mesh << yComp << flush
249 <<
"window_geometry 403 0 400 350 "
250 <<
"window_title 'Y component'" << endl;
251 z_sock <<
"solution\n" << mesh << zComp << flush
252 <<
"window_geometry 806 0 400 350 "
253 <<
"window_title 'Z component'" << endl;
261 dy_sock <<
"solution\n" << mesh << dyComp << flush
262 <<
"window_geometry 403 375 400 350 "
263 <<
"window_title 'Y component of Curl'" << endl;
264 dz_sock <<
"solution\n" << mesh << dzComp << flush
265 <<
"window_geometry 806 375 400 350 "
266 <<
"window_title 'Z component of Curl'" << endl;
276 xyMat(0,0) = 1.0; xyMat(1,1) = 1.0;
278 Vector zVec(3); zVec = 0.0; zVec(2) = 1;
303 xy_sock.precision(8);
304 xy_sock <<
"solution\n" << mesh << xyComp
305 <<
"window_title 'XY components'\n" << flush;
306 z_sock <<
"solution\n" << mesh << zComp << flush
307 <<
"window_geometry 403 0 400 350 "
308 <<
"window_title 'Z component'" << endl;
316 dxy_sock <<
"solution\n" << mesh << dxyComp << flush
317 <<
"window_geometry 0 375 400 350 "
318 <<
"window_title 'XY components of Curl'" << endl;
319 dz_sock <<
"solution\n" << mesh << dzComp << flush
320 <<
"window_geometry 403 375 400 350 "
321 <<
"window_title 'Z component of Curl'" << endl;
336 sol_sock.precision(8);
337 sol_sock <<
"solution\n" << mesh << sol
338 <<
"window_title 'Solution'" << flush << endl;
339 dsol_sock <<
"solution\n" << mesh << dsol << flush
340 <<
"window_geometry 0 375 400 350 "
341 <<
"window_title 'Curl of solution'" << endl;
356 E(0) = 1.1 *
sin(
kappa * x(0) + 0.0 * M_PI);
357 E(1) = 1.2 *
sin(
kappa * x(0) + 0.4 * M_PI);
358 E(2) = 1.3 *
sin(
kappa * x(0) + 0.9 * M_PI);
362 E(0) = 1.1 *
sin(
kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.0 * M_PI);
363 E(1) = 1.2 *
sin(
kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.4 * M_PI);
364 E(2) = 1.3 *
sin(
kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.9 * M_PI);
368 E(0) = 1.1 *
sin(
kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.0 * M_PI);
369 E(1) = 1.2 *
sin(
kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.4 * M_PI);
370 E(2) = 1.3 *
sin(
kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.9 * M_PI);
379 double c4 =
cos(
kappa * x(0) + 0.4 * M_PI);
380 double c9 =
cos(
kappa * x(0) + 0.9 * M_PI);
389 double c0 =
cos(
kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.0 * M_PI);
390 double c4 =
cos(
kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.4 * M_PI);
391 double c9 =
cos(
kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.9 * M_PI);
395 dE(2) = 1.2 * c4 - 1.1 * c0;
396 dE *=
kappa * M_SQRT1_2;
400 double s0 =
sin(
kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.0 * M_PI);
401 double c0 =
cos(
kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.0 * M_PI);
402 double s4 =
sin(
kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.4 * M_PI);
403 double c4 =
cos(
kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.4 * M_PI);
404 double c9 =
cos(
kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.9 * M_PI);
408 dE(0) = 1.2 * s4 * sk + 1.3 * M_SQRT1_2 * c9 * ck;
409 dE(1) = -1.1 * s0 * sk - 1.3 * M_SQRT1_2 * c9 * ck;
410 dE(2) = -M_SQRT1_2 * (1.1 * c0 - 1.2 * c4) * ck;
419 double s0 =
sin(
kappa * x(0) + 0.0 * M_PI);
420 double s4 =
sin(
kappa * x(0) + 0.4 * M_PI);
421 double s9 =
sin(
kappa * x(0) + 0.9 * M_PI);
423 f(0) = 2.2 * s0 + 1.2 * M_SQRT1_2 * s4;
425 M_SQRT1_2 * (1.1 * s0 + 1.3 * s9);
426 f(2) = 1.3 * (2.0 +
kappa *
kappa) * s9 + 1.2 * M_SQRT1_2 * s4;
430 double s0 =
sin(
kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.0 * M_PI);
431 double s4 =
sin(
kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.4 * M_PI);
432 double s9 =
sin(
kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.9 * M_PI);
436 f(1) = 0.55 * (M_SQRT2 - kappa *
kappa) * s0 +
437 0.6 * (4.0 + kappa * kappa) * s4 +
439 f(2) = 0.6 * M_SQRT2 * s4 + 1.3 * (2.0 + kappa *
kappa) * s9;
443 double s0 =
sin(kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.0 * M_PI);
444 double c0 =
cos(kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.0 * M_PI);
445 double s4 =
sin(kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.4 * M_PI);
446 double c4 =
cos(kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.4 * M_PI);
447 double s9 =
sin(kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.9 * M_PI);
448 double c9 =
cos(kappa * M_SQRT1_2 * (x(0) + x(1)) + 0.9 * M_PI);
449 double sk =
sin(kappa * x(2));
450 double ck =
cos(kappa * x(2));
452 f(0) = 0.55 * (4.0 + 3.0 * kappa *
kappa) * s0 * ck +
453 0.6 * (M_SQRT2 - kappa * kappa) * s4 * ck -
454 0.65 * M_SQRT2 * kappa * kappa * c9 * sk;
456 f(1) = 0.55 * (M_SQRT2 - kappa *
kappa) * s0 * ck +
457 0.6 * (4.0 + 3.0 * kappa * kappa) * s4 * ck +
458 0.65 * M_SQRT2 * s9 * ck -
459 0.65 * M_SQRT2 * kappa * kappa * c9 * sk;
461 f(2) = 0.6 * M_SQRT2 * s4 * ck -
462 M_SQRT2 * kappa * kappa * (0.55 * c0 + 0.6 * c4) * sk
463 + 1.3 * (2.0 + kappa * kappa) * s9 * ck;
A matrix coefficient that is constant in space and time.
int Size() const
Return the logical size of the array.
Class for grid function - Vector with associated FE space.
A coefficient that is constant across space and time.
Scalar coefficient defined as the inner product of two vector coefficients.
Arbitrary order 3D H(curl)-conforming Nedelec finite elements in 1D.
Integrator for (curl u, curl v) for Nedelec elements.
Vector coefficient that is constant in space and time.
Pointer to an Operator of a specified type.
Data type dense matrix using column-major storage.
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 ...
Arbitrary order 3D H(curl)-conforming Nedelec finite elements in 2D.
Data type for Gauss-Seidel smoother of sparse matrix.
Direct sparse solver using UMFPACK.
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Vector coefficient defined as the Curl of a vector GridFunction.
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)
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
FDualNumber< tbase > cos(const FDualNumber< tbase > &f)
cos([dual number])
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
A general vector function coefficient.
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
double Control[UMFPACK_CONTROL]
FDualNumber< tbase > sin(const FDualNumber< tbase > &f)
sin([dual number])
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
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 f_exact(const Vector &, Vector &)
Vector coefficient defined as a product of a matrix coefficient and a vector coefficient.
virtual double ComputeHCurlError(VectorCoefficient *exsol, VectorCoefficient *excurl, const IntegrationRule *irs[]=NULL) const
Returns the error measured in H(curl)-norm for ND elements.
virtual void Print(std::ostream &os=mfem::out) const
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
void E_exact(const Vector &, Vector &)
for VectorFiniteElements (Nedelec, Raviart-Thomas)
Arbitrary order H(curl)-conforming Nedelec finite elements.
Vector coefficient defined by a vector GridFunction.
void CurlE_exact(const Vector &, Vector &)
Arbitrary order H1-conforming (continuous) finite elements.
void ParseCheck(std::ostream &out=mfem::out)
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Arbitrary order "L2-conforming" discontinuous finite elements.
double f(const Vector &p)
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a SparseMatrix.
double sigma(const Vector &x)