MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
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
28using namespace std;
29using namespace mfem;
30
31Mesh * GetMesh(int type);
32
33void trans(const Vector &x, Vector &r);
34
35void sigmaFunc(const Vector &x, DenseMatrix &s);
36
38{
39 return (0.25 * (2.0 + x[0]) - x[2]) * (x[2] + 0.25 * (2.0 + x[0]));
40}
41
42void 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
50void 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
64int 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);
130 b.AddDomainIntegrator(new DomainLFIntegrator(one));
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);
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 real_t 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 real_t 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.
218Mesh * 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.
303void trans(const Vector &x, Vector &r)
304{
305 r.SetSize(3);
306
307 real_t tol = 1e-6;
308 real_t 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
338{
339 s.SetSize(3);
340 real_t 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}
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
int Size() const
Return the logical size of the array.
Definition array.hpp:144
Abstract base class BilinearFormIntegrator.
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
A coefficient that is constant across space and time.
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
Class for domain integration .
Definition lininteg.hpp:109
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition fespace.hpp:716
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1) const
Get a list of essential true dofs, ess_tdof_list, corresponding to the boundary attributes marked in ...
Definition fespace.cpp:588
A general function coefficient.
Data type for Gauss-Seidel smoother of sparse matrix.
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
virtual real_t ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const
virtual void ComputeFlux(BilinearFormIntegrator &blfi, GridFunction &flux, bool wcoef=true, int subdomain=-1)
Definition gridfunc.cpp:308
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
Vector with associated FE space and LinearFormIntegrators.
A matrix coefficient with an optional scalar coefficient multiplier q. The matrix function can either...
Mesh data type.
Definition mesh.hpp:56
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition mesh.hpp:282
int AddQuad(int v1, int v2, int v3, int v4, int attr=1)
Adds a quadrilateral to the mesh given by 4 vertices v1 through v4.
Definition mesh.cpp:1743
int AddTriangle(int v1, int v2, int v3, int attr=1)
Adds a triangle to the mesh given by 3 vertices v1 through v3.
Definition mesh.cpp:1729
void FinalizeTopology(bool generate_bdr=true)
Finalize the construction of the secondary topology (connectivity) data of a Mesh.
Definition mesh.cpp:3135
virtual void Print(std::ostream &os=mfem::out, const std::string &comments="") const
Definition mesh.hpp:2288
int AddVertex(real_t x, real_t y=0.0, real_t z=0.0)
Definition mesh.cpp:1658
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
int AddBdrSegment(int v1, int v2, int attr=1)
Definition mesh.cpp:2035
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:6211
void Transform(void(*f)(const Vector &, Vector &))
Definition mesh.cpp:12821
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
Pointer to an Operator of a specified type.
Definition handle.hpp:34
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition operator.hpp:66
void ParseCheck(std::ostream &out=mfem::out)
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...
Definition optparser.hpp:82
Data type sparse matrix.
Definition sparsemat.hpp:51
A general vector function coefficient.
Vector data type.
Definition vector.hpp:80
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
real_t sigma(const Vector &x)
Definition maxwell.cpp:102
int dim
Definition ex24.cpp:53
void sigmaFunc(const Vector &x, DenseMatrix &s)
Definition ex29.cpp:337
void duExact(const Vector &x, Vector &du)
Definition ex29.cpp:42
void fluxExact(const Vector &x, Vector &f)
Definition ex29.cpp:50
real_t uExact(const Vector &x)
Definition ex29.cpp:37
Mesh * GetMesh(int type)
Definition ex29.cpp:218
void trans(const Vector &x, Vector &r)
Definition ex29.cpp:303
int main()
real_t b
Definition lissajous.cpp:42
real_t a
Definition lissajous.cpp:41
const int visport
void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x, int print_iter, int max_num_iter, real_t RTOLERANCE, real_t ATOLERANCE)
Preconditioned conjugate gradient method. (tolerances are squared)
Definition solvers.cpp:913
float real_t
Definition config.hpp:43
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30
const char vishost[]
RefCoord s[3]