MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
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);
75  args.AddOption(&mesh_type, "-mt", "--mesh-type",
76  "Mesh type: 3 - Triangular, 4 - Quadrilateral.");
77  args.AddOption(&mesh_order, "-mo", "--mesh-order",
78  "Geometric order of the curved mesh.");
79  args.AddOption(&ref_levels, "-r", "--refine",
80  "Number of times to refine the mesh uniformly in serial.");
81  args.AddOption(&order, "-o", "--order",
82  "Finite element order (polynomial degree).");
83  args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
84  "--no-static-condensation", "Enable static condensation.");
85  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
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);
144  BilinearFormIntegrator *integ = new DiffusionIntegrator(sigma);
145  a.AddDomainIntegrator(integ);
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 
226  mesh->AddVertex(-1.0, -1.0, 0.0);
227  mesh->AddVertex( 1.0, -1.0, 0.0);
228  mesh->AddVertex( 1.0, 1.0, 0.0);
229  mesh->AddVertex(-1.0, 1.0, 0.0);
230  mesh->AddVertex(-1.0, -1.0, 1.0);
231  mesh->AddVertex( 1.0, -1.0, 1.0);
232  mesh->AddVertex( 1.0, 1.0, 1.0);
233  mesh->AddVertex(-1.0, 1.0, 1.0);
234  mesh->AddVertex( 0.0, -1.0, 0.5);
235  mesh->AddVertex( 1.0, 0.0, 0.5);
236  mesh->AddVertex( 0.0, 1.0, 0.5);
237  mesh->AddVertex(-1.0, 0.0, 0.5);
238 
239  mesh->AddTriangle(0, 1, 8);
240  mesh->AddTriangle(1, 5, 8);
241  mesh->AddTriangle(5, 4, 8);
242  mesh->AddTriangle(4, 0, 8);
243  mesh->AddTriangle(1, 2, 9);
244  mesh->AddTriangle(2, 6, 9);
245  mesh->AddTriangle(6, 5, 9);
246  mesh->AddTriangle(5, 1, 9);
247  mesh->AddTriangle(2, 3, 10);
248  mesh->AddTriangle(3, 7, 10);
249  mesh->AddTriangle(7, 6, 10);
250  mesh->AddTriangle(6, 2, 10);
251  mesh->AddTriangle(3, 0, 11);
252  mesh->AddTriangle(0, 4, 11);
253  mesh->AddTriangle(4, 7, 11);
254  mesh->AddTriangle(7, 3, 11);
255 
256  mesh->AddBdrSegment(0, 1, 1);
257  mesh->AddBdrSegment(1, 2, 1);
258  mesh->AddBdrSegment(2, 3, 1);
259  mesh->AddBdrSegment(3, 0, 1);
260  mesh->AddBdrSegment(5, 4, 2);
261  mesh->AddBdrSegment(6, 5, 2);
262  mesh->AddBdrSegment(7, 6, 2);
263  mesh->AddBdrSegment(4, 7, 2);
264  }
265  else if (type == 4)
266  {
267  mesh = new Mesh(2, 8, 4, 8, 3);
268 
269  mesh->AddVertex(-1.0, -1.0, 0.0);
270  mesh->AddVertex( 1.0, -1.0, 0.0);
271  mesh->AddVertex( 1.0, 1.0, 0.0);
272  mesh->AddVertex(-1.0, 1.0, 0.0);
273  mesh->AddVertex(-1.0, -1.0, 1.0);
274  mesh->AddVertex( 1.0, -1.0, 1.0);
275  mesh->AddVertex( 1.0, 1.0, 1.0);
276  mesh->AddVertex(-1.0, 1.0, 1.0);
277 
278  mesh->AddQuad(0, 1, 5, 4);
279  mesh->AddQuad(1, 2, 6, 5);
280  mesh->AddQuad(2, 3, 7, 6);
281  mesh->AddQuad(3, 0, 4, 7);
282 
283  mesh->AddBdrSegment(0, 1, 1);
284  mesh->AddBdrSegment(1, 2, 1);
285  mesh->AddBdrSegment(2, 3, 1);
286  mesh->AddBdrSegment(3, 0, 1);
287  mesh->AddBdrSegment(5, 4, 2);
288  mesh->AddBdrSegment(6, 5, 2);
289  mesh->AddBdrSegment(7, 6, 2);
290  mesh->AddBdrSegment(4, 7, 2);
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 }
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:108
void duExact(const Vector &x, Vector &du)
Definition: ex29.cpp:42
void trans(const Vector &u, Vector &x)
Definition: ex27.cpp:412
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
int AddQuad(int v1, int v2, int v3, int v4, int attr=1)
Definition: mesh.cpp:1674
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:513
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:168
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:1660
virtual void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
Form the linear system A X = B, corresponding to this bilinear form and the linear form b(...
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
Recover the solution of a linear system formed with FormLinearSystem().
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
void Transform(void(*f)(const Vector &, Vector &))
Definition: mesh.cpp:11629
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:564
Data type for Gauss-Seidel smoother of sparse matrix.
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:3676
int AddVertex(double x, double y=0.0, double z=0.0)
Definition: mesh.cpp:1613
Data type sparse matrix.
Definition: sparsemat.hpp:50
A matrix coefficient with an optional scalar coefficient multiplier q. The matrix function can either...
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:67
constexpr char vishost[]
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9816
constexpr int visport
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:5343
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
int AddBdrSegment(int v1, int v2, int attr=1)
Definition: mesh.cpp:1823
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:904
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:590
int Dimension() const
Definition: mesh.hpp:1006
A general vector function coefficient.
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:41
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:35
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:270
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
void FinalizeTopology(bool generate_bdr=true)
Finalize the construction of the secondary topology (connectivity) data of a Mesh.
Definition: mesh.cpp:2890
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
virtual double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const
Definition: gridfunc.cpp:2766
double a
Definition: lissajous.cpp:41
A &quot;square matrix&quot; operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
virtual void Print(std::ostream &os=mfem::out) const
Definition: mesh.hpp:1671
double uExact(const Vector &x)
Definition: ex29.cpp:37
Mesh * GetMesh(int type)
Definition: ex29.cpp:218
int dim
Definition: ex24.cpp:53
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
A general function coefficient.
Vector data type.
Definition: vector.hpp:60
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:173
virtual void ComputeFlux(BilinearFormIntegrator &blfi, GridFunction &flux, bool wcoef=true, int subdomain=-1)
Definition: gridfunc.cpp:309
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:220
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
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
int main()
void ParseCheck(std::ostream &out=mfem::out)
Definition: optparser.cpp:252
void EnableStaticCondensation()
Enable the use of static condensation. For details see the description for class StaticCondensation i...
double f(const Vector &p)
double sigma(const Vector &x)
Definition: maxwell.cpp:102