MFEM  v4.6.0 Finite element discretization library
ex29.cpp
Go to the documentation of this file.
1 // MFEM Example 29
2 //
3 // Compile with: make ex29
4 //
5 // Sample runs: ex29
6 // ex29 -r 2 -sc
7 // ex29 -mt 3 -o 4 -sc
8 // ex29 -mt 3 -r 2 -o 4 -sc
9 //
10 // Description: This example code demonstrates the use of MFEM to define a
11 // finite element discretization of a PDE on a 2 dimensional
12 // surface embedded in a 3 dimensional domain. In this case we
13 // solve the Laplace problem -Div(sigma Grad u) = 1, with
14 // homogeneous Dirichlet boundary conditions, where sigma is an
15 // anisotropic diffusion constant defined as a 3x3 matrix
16 // coefficient.
17 //
18 // This example demonstrates the use of finite element integrators
19 // on 2D domains with 3D coefficients.
20 //
21 // We recommend viewing examples 1 and 7 before viewing this
22 // example.
23
24 #include "mfem.hpp"
25 #include <fstream>
26 #include <iostream>
27
28 using namespace std;
29 using namespace mfem;
30
31 Mesh * GetMesh(int type);
32
33 void trans(const Vector &x, Vector &r);
34
35 void sigmaFunc(const Vector &x, DenseMatrix &s);
36
37 double uExact(const Vector &x)
38 {
39  return (0.25 * (2.0 + x[0]) - x[2]) * (x[2] + 0.25 * (2.0 + x[0]));
40 }
41
42 void duExact(const Vector &x, Vector &du)
43 {
44  du.SetSize(3);
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];
47  du[2] = -2.0 * x[2];
48 }
49
50 void fluxExact(const Vector &x, Vector &f)
51 {
52  f.SetSize(3);
53
54  DenseMatrix s(3);
55  sigmaFunc(x, s);
56
57  Vector du(3);
58  duExact(x, du);
59
60  s.Mult(du, f);
61  f *= -1.0;
62 }
63
64 int main(int argc, char *argv[])
65 {
66  // 1. Parse command-line options.
67  int order = 3;
68  int mesh_type = 4; // Default to Quadrilateral mesh
69  int mesh_order = 3;
70  int ref_levels = 0;
71  bool static_cond = false;
72  bool visualization = true;
73
74  OptionsParser args(argc, argv);
76  "Mesh type: 3 - Triangular, 4 - Quadrilateral.");
78  "Geometric order of the curved mesh.");
80  "Number of times to refine the mesh uniformly in serial.");
82  "Finite element order (polynomial degree).");
84  "--no-static-condensation", "Enable static condensation.");
86  "--no-visualization",
87  "Enable or disable GLVis visualization.");
88  args.ParseCheck();
89
90  // 2. Construct a quadrilateral or triangular mesh with the topology of a
91  // cylindrical surface.
92  Mesh *mesh = GetMesh(mesh_type);
93  int dim = mesh->Dimension();
94
95  // 3. Refine the mesh to increase the resolution. In this example we do
96  // 'ref_levels' of uniform refinement.
97  for (int l = 0; l < ref_levels; l++)
98  {
99  mesh->UniformRefinement();
100  }
101
102  // 4. Transform the mesh so that it has a more interesting geometry.
103  mesh->SetCurvature(mesh_order);
104  mesh->Transform(trans);
105
106  // 5. Define a finite element space on the mesh. Here we use continuous
107  // Lagrange finite elements of the specified order.
108  H1_FECollection fec(order, dim);
109  FiniteElementSpace fespace(mesh, &fec);
110  cout << "Number of finite element unknowns: "
111  << fespace.GetTrueVSize() << endl;
112
113  // 6. Determine the list of true (i.e. conforming) essential boundary dofs.
114  // In this example, the boundary conditions are defined by marking all
115  // the boundary attributes from the mesh as essential (Dirichlet) and
116  // converting them to a list of true dofs.
117  Array<int> ess_tdof_list;
118  if (mesh->bdr_attributes.Size())
119  {
120  Array<int> ess_bdr(mesh->bdr_attributes.Max());
121  ess_bdr = 1;
122  fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
123  }
124
125  // 7. Set up the linear form b(.) which corresponds to the right-hand side of
126  // the FEM linear system, which in this case is (1,phi_i) where phi_i are
127  // the basis functions in the finite element fespace.
128  LinearForm b(&fespace);
129  ConstantCoefficient one(1.0);
131  b.Assemble();
132
133  // 8. Define the solution vector x as a finite element grid function
134  // corresponding to fespace. Initialize x with initial guess of zero,
135  // which satisfies the boundary conditions.
136  GridFunction x(&fespace);
137  x = 0.0;
138
139  // 9. Set up the bilinear form a(.,.) on the finite element space
140  // corresponding to the Laplacian operator -Delta, by adding the Diffusion
141  // domain integrator.
142  BilinearForm a(&fespace);
146
147  // 10. Assemble the bilinear form and the corresponding linear system,
148  // applying any necessary transformations such as: eliminating boundary
149  // conditions, applying conforming constraints for non-conforming AMR,
150  // static condensation, etc.
151  if (static_cond) { a.EnableStaticCondensation(); }
152  a.Assemble();
153
154  OperatorPtr A;
155  Vector B, X;
156  a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
157
158  cout << "Size of linear system: " << A->Height() << endl;
159
160  // 11. Solve the linear system A X = B.
161  // Use a simple symmetric Gauss-Seidel preconditioner with PCG.
162  GSSmoother M((SparseMatrix&)(*A));
163  PCG(*A, M, B, X, 1, 200, 1e-12, 0.0);
164
165  // 12. Recover the solution as a finite element grid function.
166  a.RecoverFEMSolution(X, b, x);
167
168  // 13. Compute error in the solution and its flux
170  double error = x.ComputeL2Error(uCoef);
171
172  cout << "|u - u_h|_2 = " << error << endl;
173
174  FiniteElementSpace flux_fespace(mesh, &fec, 3);
175  GridFunction flux(&flux_fespace);
176  x.ComputeFlux(*integ, flux); flux *= -1.0;
177
179  double flux_err = flux.ComputeL2Error(fluxCoef);
180
181  cout << "|f - f_h|_2 = " << flux_err << endl;
182
183  // 14. Save the refined mesh and the solution. This output can be viewed
184  // later using GLVis: "glvis -m refined.mesh -g sol.gf".
185  ofstream mesh_ofs("refined.mesh");
186  mesh_ofs.precision(8);
187  mesh->Print(mesh_ofs);
188  ofstream sol_ofs("sol.gf");
189  sol_ofs.precision(8);
190  x.Save(sol_ofs);
191
192  // 15. Send the solution by socket to a GLVis server.
193  if (visualization)
194  {
195  char vishost[] = "localhost";
196  int visport = 19916;
197  socketstream sol_sock(vishost, visport);
198  sol_sock.precision(8);
199  sol_sock << "solution\n" << *mesh << x
200  << "window_title 'Solution'\n" << flush;
201
202  socketstream flux_sock(vishost, visport);
203  flux_sock.precision(8);
204  flux_sock << "solution\n" << *mesh << flux
205  << "keys vvv\n"
206  << "window_geometry 402 0 400 350\n"
207  << "window_title 'Flux'\n" << flush;
208  }
209
210  // 16. Free the used memory.
211  delete mesh;
212
213  return 0;
214 }
215
216 // Defines a mesh consisting of four flat rectangular surfaces connected to form
217 // a loop.
218 Mesh * GetMesh(int type)
219 {
220  Mesh * mesh = NULL;
221
222  if (type == 3)
223  {
224  mesh = new Mesh(2, 12, 16, 8, 3);
225
238
255
264  }
265  else if (type == 4)
266  {
267  mesh = new Mesh(2, 8, 4, 8, 3);
268
277
282
291  }
292  else
293  {
294  MFEM_ABORT("Unrecognized mesh type " << type << "!");
295  }
296  mesh->FinalizeTopology();
297
298  return mesh;
299 }
300
301 // Transforms the four-sided loop into a curved cylinder with skewed top and
302 // base.
303 void trans(const Vector &x, Vector &r)
304 {
305  r.SetSize(3);
306
307  double tol = 1e-6;
308  double theta = 0.0;
309  if (fabs(x[1] + 1.0) < tol)
310  {
311  theta = 0.25 * M_PI * (x[0] - 2.0);
312  }
313  else if (fabs(x[0] - 1.0) < tol)
314  {
315  theta = 0.25 * M_PI * x[1];
316  }
317  else if (fabs(x[1] - 1.0) < tol)
318  {
319  theta = 0.25 * M_PI * (2.0 - x[0]);
320  }
321  else if (fabs(x[0] + 1.0) < tol)
322  {
323  theta = 0.25 * M_PI * (4.0 - x[1]);
324  }
325  else
326  {
327  cerr << "side not recognized "
328  << x[0] << " " << x[1] << " " << x[2] << endl;
329  }
330
331  r[0] = cos(theta);
332  r[1] = sin(theta);
333  r[2] = 0.25 * (2.0 * x[2] - 1.0) * (r[0] + 2.0);
334 }
335
336 // Anisotropic diffusion coefficient
337 void sigmaFunc(const Vector &x, DenseMatrix &s)
338 {
339  s.SetSize(3);
340  double a = 17.0 - 2.0 * x[0] * (1.0 + x[0]);
341  s(0,0) = 0.5 + x[0] * x[0] * (8.0 / a - 0.5);
342  s(0,1) = x[0] * x[1] * (8.0 / a - 0.5);
343  s(0,2) = 0.0;
344  s(1,0) = s(0,1);
345  s(1,1) = 0.5 * x[0] * x[0] + 8.0 * x[1] * x[1] / a;
346  s(1,2) = 0.0;
347  s(2,0) = 0.0;
348  s(2,1) = 0.0;
349  s(2,2) = a / 32.0;
350 }
virtual double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const
Definition: gridfunc.cpp:2786
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:108
void duExact(const Vector &x, Vector &du)
Definition: ex29.cpp:42
int visport
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
Definition: mesh.cpp:1694
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:517
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
void sigmaFunc(const Vector &x, DenseMatrix &s)
Definition: ex29.cpp:337
int AddTriangle(int v1, int v2, int v3, int attr=1)
Definition: mesh.cpp:1680
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:68
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
void Transform(void(*f)(const Vector &, Vector &))
Definition: mesh.cpp:12045
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 ...
Definition: fespace.cpp:587
STL namespace.
Data type for Gauss-Seidel smoother of sparse matrix.
int AddVertex(double x, double y=0.0, double z=0.0)
Definition: mesh.cpp:1626
char vishost[]
Data type sparse matrix.
Definition: sparsemat.hpp:50
A matrix coefficient with an optional scalar coefficient multiplier q. The matrix function can either...
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:10232
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Set the curvature of the mesh nodes using the given polynomial degree.
Definition: mesh.cpp:5635
int AddBdrSegment(int v1, int v2, int attr=1)
Definition: mesh.cpp:1843
int main(int argc, char *argv[])
Definition: ex29.cpp:64
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)
Definition: solvers.cpp:913
A general vector function coefficient.
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:712
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:24
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:275
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:219
void FinalizeTopology(bool generate_bdr=true)
Finalize the construction of the secondary topology (connectivity) data of a Mesh.
Definition: mesh.cpp:2945
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 &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
void trans(const Vector &x, Vector &r)
Definition: ex29.cpp:303
double a
Definition: lissajous.cpp:41
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
double uExact(const Vector &x)
Definition: ex29.cpp:37
Mesh * GetMesh(int type)
Definition: ex29.cpp:218
int dim
Definition: ex24.cpp:53
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
A general function coefficient.
Vector data type.
Definition: vector.hpp:58
virtual void Print(std::ostream &os=mfem::out) const
Definition: mesh.hpp:2011
virtual void ComputeFlux(BilinearFormIntegrator &blfi, GridFunction &flux, bool wcoef=true, int subdomain=-1)
Definition: gridfunc.cpp:308
double sigma(const Vector &x)
Definition: maxwell.cpp:102
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:259
Vector with associated FE space and LinearFormIntegrators.
Definition: linearform.hpp:24
RefCoord s[3]
void fluxExact(const Vector &x, Vector &f)
Definition: ex29.cpp:50
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:3696
void ParseCheck(std::ostream &out=mfem::out)
Definition: optparser.cpp:255
double f(const Vector &p)