MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
ex17p.cpp
Go to the documentation of this file.
1// MFEM Example 17 - Parallel Version
2//
3// Compile with: make ex17p
4//
5// Sample runs:
6//
7// mpirun -np 4 ex17p -m ../data/beam-tri.mesh
8// mpirun -np 4 ex17p -m ../data/beam-quad.mesh
9// mpirun -np 4 ex17p -m ../data/beam-tet.mesh
10// mpirun -np 4 ex17p -m ../data/beam-hex.mesh
11// mpirun -np 4 ex17p -m ../data/beam-wedge.mesh
12// mpirun -np 4 ex17p -m ../data/beam-quad.mesh -rs 2 -rp 2 -o 3 -elast
13// mpirun -np 4 ex17p -m ../data/beam-quad.mesh -rs 2 -rp 3 -o 2 -a 1 -k 1
14// mpirun -np 4 ex17p -m ../data/beam-hex.mesh -rs 2 -rp 1 -o 2
15//
16// Description: This example code solves a simple linear elasticity problem
17// describing a multi-material cantilever beam using symmetric or
18// non-symmetric discontinuous Galerkin (DG) formulation.
19//
20// Specifically, we approximate the weak form of -div(sigma(u))=0
21// where sigma(u)=lambda*div(u)*I+mu*(grad*u+u*grad) is the stress
22// tensor corresponding to displacement field u, and lambda and mu
23// are the material Lame constants. The boundary conditions are
24// Dirichlet, u=u_D on the fixed part of the boundary, namely
25// boundary attributes 1 and 2; on the rest of the boundary we use
26// sigma(u).n=0 b.c. The geometry of the domain is assumed to be
27// as follows:
28//
29// +----------+----------+
30// boundary --->| material | material |<--- boundary
31// attribute 1 | 1 | 2 | attribute 2
32// (fixed) +----------+----------+ (fixed, nonzero)
33//
34// The example demonstrates the use of high-order DG vector finite
35// element spaces with the linear DG elasticity bilinear form,
36// meshes with curved elements, and the definition of piece-wise
37// constant and function vector-coefficient objects. The use of
38// non-homogeneous Dirichlet b.c. imposed weakly, is also
39// illustrated.
40//
41// We recommend viewing examples 2p and 14p before viewing this
42// example.
43
44#include "mfem.hpp"
45#include <fstream>
46#include <iostream>
47
48using namespace std;
49using namespace mfem;
50
51// Initial displacement, used for Dirichlet boundary conditions on boundary
52// attributes 1 and 2.
53void InitDisplacement(const Vector &x, Vector &u);
54
55// A Coefficient for computing the components of the stress.
56class StressCoefficient : public Coefficient
57{
58protected:
59 Coefficient &lambda, &mu;
60 GridFunction *u; // displacement
61 int si, sj; // component of the stress to evaluate, 0 <= si,sj < dim
62
63 DenseMatrix grad; // auxiliary matrix, used in Eval
64
65public:
66 StressCoefficient(Coefficient &lambda_, Coefficient &mu_)
67 : lambda(lambda_), mu(mu_), u(NULL), si(0), sj(0) { }
68
69 void SetDisplacement(GridFunction &u_) { u = &u_; }
70 void SetComponent(int i, int j) { si = i; sj = j; }
71
72 virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip);
73};
74
75// Simple GLVis visualization manager.
76class VisMan : public iostream
77{
78protected:
79 const char *host;
80 int port;
82 int sid; // active socket, index inside 'sock'.
83
84 int win_x, win_y, win_w, win_h;
85 int win_stride_x, win_stride_y, win_nx;
86
87public:
88 VisMan(const char *vishost, const int visport);
89 void NewWindow();
90 void CloseConnection();
91 void PositionWindow();
92 virtual ~VisMan();
93};
94
95// Manipulators for the GLVis visualization manager.
96void new_window (VisMan &v) { v.NewWindow(); }
97void position_window (VisMan &v) { v.PositionWindow(); }
98void close_connection(VisMan &v) { v.CloseConnection(); }
99ostream &operator<<(ostream &v, void (*f)(VisMan&));
100
101int main(int argc, char *argv[])
102{
103 Mpi::Init(argc, argv);
104 Hypre::Init();
105
106 // 1. Define and parse command-line options.
107 const char *mesh_file = "../data/beam-tri.mesh";
108 int ser_ref_levels = -1;
109 int par_ref_levels = 1;
110 int order = 1;
111 real_t alpha = -1.0;
112 real_t kappa = -1.0;
113 bool amg_elast = false;
114 bool visualization = 1;
115
116 OptionsParser args(argc, argv);
117 args.AddOption(&mesh_file, "-m", "--mesh",
118 "Mesh file to use.");
119 args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
120 "Number of times to refine the mesh uniformly before parallel"
121 " partitioning, -1 for auto.");
122 args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
123 "Number of times to refine the mesh uniformly after parallel"
124 " partitioning.");
125 args.AddOption(&order, "-o", "--order",
126 "Finite element order (polynomial degree).");
127 args.AddOption(&alpha, "-a", "--alpha",
128 "One of the two DG penalty parameters, typically +1/-1."
129 " See the documentation of class DGElasticityIntegrator.");
130 args.AddOption(&kappa, "-k", "--kappa",
131 "One of the two DG penalty parameters, should be positive."
132 " Negative values are replaced with (order+1)^2.");
133 args.AddOption(&amg_elast, "-elast", "--amg-for-elasticity", "-sys",
134 "--amg-for-systems",
135 "Use the special AMG elasticity solver (GM/LN approaches), "
136 "or standard AMG for systems (unknown approach).");
137 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
138 "--no-visualization",
139 "Enable or disable GLVis visualization.");
140 args.Parse();
141 if (!args.Good())
142 {
143 if (Mpi::Root()) { args.PrintUsage(cout); }
144 return 1;
145 }
146 if (kappa < 0)
147 {
148 kappa = (order+1)*(order+1);
149 }
150 if (Mpi::Root()) { args.PrintOptions(cout); }
151
152 // 2. Read the mesh from the given mesh file.
153 Mesh mesh(mesh_file, 1, 1);
154 int dim = mesh.Dimension();
155
156 if (mesh.attributes.Max() < 2 || mesh.bdr_attributes.Max() < 2)
157 {
158 if (Mpi::Root())
159 {
160 cerr << "\nInput mesh should have at least two materials and "
161 << "two boundary attributes! (See schematic in ex17p.cpp)\n"
162 << endl;
163 }
164 return 3;
165 }
166
167 // 3. Refine the mesh to increase the resolution.
168 if (ser_ref_levels < 0)
169 {
170 ser_ref_levels = (int)floor(log(5000./mesh.GetNE())/log(2.)/dim);
171 }
172 for (int l = 0; l < ser_ref_levels; l++)
173 {
174 mesh.UniformRefinement();
175 }
176 // Since NURBS meshes do not support DG integrators, we convert them to
177 // regular polynomial mesh of the specified (solution) order.
178 if (mesh.NURBSext) { mesh.SetCurvature(order); }
179
180 ParMesh pmesh(MPI_COMM_WORLD, mesh);
181 mesh.Clear();
182
183 for (int l = 0; l < par_ref_levels; l++)
184 {
185 pmesh.UniformRefinement();
186 }
187
188 // 4. Define a DG vector finite element space on the mesh. Here, we use
189 // Gauss-Lobatto nodal basis because it gives rise to a sparser matrix
190 // compared to the default Gauss-Legendre nodal basis.
192 ParFiniteElementSpace fespace(&pmesh, &fec, dim, Ordering::byVDIM);
193
194 HYPRE_BigInt glob_size = fespace.GlobalTrueVSize();
195 if (Mpi::Root())
196 {
197 cout << "Number of finite element unknowns: " << glob_size
198 << "\nAssembling: " << flush;
199 }
200
201 // 5. In this example, the Dirichlet boundary conditions are defined by
202 // marking boundary attributes 1 and 2 in the marker Array 'dir_bdr'.
203 // These b.c. are imposed weakly, by adding the appropriate boundary
204 // integrators over the marked 'dir_bdr' to the bilinear and linear forms.
205 // With this DG formulation, there are no essential boundary conditions.
206 Array<int> ess_tdof_list; // no essential b.c. (empty list)
207 Array<int> dir_bdr(pmesh.bdr_attributes.Max());
208 dir_bdr = 0;
209 dir_bdr[0] = 1; // boundary attribute 1 is Dirichlet
210 dir_bdr[1] = 1; // boundary attribute 2 is Dirichlet
211
212 // 6. Define the DG solution vector 'x' as a finite element grid function
213 // corresponding to fespace. Initialize 'x' using the 'InitDisplacement'
214 // function.
215 ParGridFunction x(&fespace);
217 x.ProjectCoefficient(init_x);
218
219 // 7. Set up the Lame constants for the two materials. They are defined as
220 // piece-wise (with respect to the element attributes) constant
221 // coefficients, i.e. type PWConstCoefficient.
222 Vector lambda(pmesh.attributes.Max());
223 lambda = 1.0; // Set lambda = 1 for all element attributes.
224 lambda(0) = 50.0; // Set lambda = 50 for element attribute 1.
225 PWConstCoefficient lambda_c(lambda);
226 Vector mu(pmesh.attributes.Max());
227 mu = 1.0; // Set mu = 1 for all element attributes.
228 mu(0) = 50.0; // Set mu = 50 for element attribute 1.
230
231 // 8. Set up the linear form b(.) which corresponds to the right-hand side of
232 // the FEM linear system. In this example, the linear form b(.) consists
233 // only of the terms responsible for imposing weakly the Dirichlet
234 // boundary conditions, over the attributes marked in 'dir_bdr'. The
235 // values for the Dirichlet boundary condition are taken from the
236 // VectorFunctionCoefficient 'x_init' which in turn is based on the
237 // function 'InitDisplacement'.
238 ParLinearForm b(&fespace);
239 if (Mpi::Root()) { cout << "r.h.s. ... " << flush; }
240 b.AddBdrFaceIntegrator(
242 init_x, lambda_c, mu_c, alpha, kappa), dir_bdr);
243 b.Assemble();
244
245 // 9. Set up the bilinear form a(.,.) on the DG finite element space
246 // corresponding to the linear elasticity integrator with coefficients
247 // lambda and mu as defined above. The additional interior face integrator
248 // ensures the weak continuity of the displacement field. The additional
249 // boundary face integrator works together with the boundary integrator
250 // added to the linear form b(.) to impose weakly the Dirichlet boundary
251 // conditions.
252 ParBilinearForm a(&fespace);
253 a.AddDomainIntegrator(new ElasticityIntegrator(lambda_c, mu_c));
254 a.AddInteriorFaceIntegrator(
255 new DGElasticityIntegrator(lambda_c, mu_c, alpha, kappa));
256 a.AddBdrFaceIntegrator(
257 new DGElasticityIntegrator(lambda_c, mu_c, alpha, kappa), dir_bdr);
258
259 // 10. Assemble the bilinear form and the corresponding linear system.
260 if (Mpi::Root()) { cout << "matrix ... " << flush; }
261 a.Assemble();
262
264 Vector B, X;
265 a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
266 if (Mpi::Root()) { cout << "done." << endl; }
267
268 // 11. Define a simple symmetric Gauss-Seidel preconditioner and use it to
269 // solve the system Ax=b with PCG for the symmetric formulation, or GMRES
270 // for the non-symmetric.
271 const real_t rtol = 1e-6;
272 HypreBoomerAMG amg(A);
273 if (amg_elast)
274 {
275 amg.SetElasticityOptions(&fespace);
276 }
277 else
278 {
280 }
281 CGSolver pcg(A.GetComm());
282 GMRESSolver gmres(A.GetComm());
283 gmres.SetKDim(50);
284 IterativeSolver &ipcg = pcg, &igmres = gmres;
285 IterativeSolver &solver = (alpha == -1.0) ? ipcg : igmres;
286 solver.SetRelTol(rtol);
287 solver.SetMaxIter(500);
288 solver.SetPrintLevel(1);
289 solver.SetOperator(A);
290 solver.SetPreconditioner(amg);
291 solver.Mult(B, X);
292
293 // 12. Recover the solution as a finite element grid function 'x'.
294 a.RecoverFEMSolution(X, b, x);
295
296 // 13. Use the DG solution space as the mesh nodal space. This allows us to
297 // save the displaced mesh as a curved DG mesh.
298 pmesh.SetNodalFESpace(&fespace);
299
300 Vector reference_nodes;
301 if (visualization) { reference_nodes = *pmesh.GetNodes(); }
302
303 // 14. Save the displaced mesh and minus the solution (which gives the
304 // backward displacements to the reference mesh). This output can be
305 // viewed later using GLVis: "glvis -m displaced.mesh -g sol.gf".
306 {
307 *pmesh.GetNodes() += x;
308 x.Neg(); // x = -x
309
310 ostringstream mesh_name, sol_name;
311 mesh_name << "mesh." << setfill('0') << setw(6) << Mpi::WorldRank();
312 sol_name << "sol." << setfill('0') << setw(6) << Mpi::WorldRank();
313
314 ofstream mesh_ofs(mesh_name.str().c_str());
315 mesh_ofs.precision(8);
316 mesh_ofs << pmesh;
317
318 ofstream sol_ofs(sol_name.str().c_str());
319 sol_ofs.precision(8);
320 sol_ofs << x;
321 }
322
323 // 15. Visualization: send data by socket to a GLVis server.
324 if (visualization)
325 {
326 char vishost[] = "localhost";
327 int visport = 19916;
328 VisMan vis(vishost, visport);
329 const char *glvis_keys = (dim < 3) ? "Rjlc" : "c";
330
331 // Visualize the deformed configuration.
332 vis << new_window << setprecision(8)
333 << "parallel " << pmesh.GetNRanks() << ' ' << pmesh.GetMyRank()
334 << '\n'
335 << "solution\n" << pmesh << x << flush
336 << "keys " << glvis_keys << endl
337 << "window_title 'Deformed configuration'" << endl
338 << "plot_caption 'Backward displacement'" << endl
340
341 // Visualize the stress components.
342 const char *c = "xyz";
343 ParFiniteElementSpace scalar_dg_space(&pmesh, &fec);
344 ParGridFunction stress(&scalar_dg_space);
345 StressCoefficient stress_c(lambda_c, mu_c);
346 *pmesh.GetNodes() = reference_nodes;
347 x.Neg(); // x = -x
348 stress_c.SetDisplacement(x);
349 for (int si = 0; si < dim; si++)
350 {
351 for (int sj = si; sj < dim; sj++)
352 {
353 stress_c.SetComponent(si, sj);
354 stress.ProjectCoefficient(stress_c);
355
356 MPI_Barrier(MPI_COMM_WORLD);
357 vis << new_window << setprecision(8)
358 << "parallel " << pmesh.GetNRanks() << ' ' << pmesh.GetMyRank()
359 << '\n'
360 << "solution\n" << pmesh << stress << flush
361 << "keys " << glvis_keys << endl
362 << "window_title |Stress " << c[si] << c[sj] << '|' << endl
364 }
365 }
366 }
367
368 return 0;
369}
370
371
373{
374 u = 0.0;
375 u(u.Size()-1) = -0.2*x(0);
376}
377
378
379real_t StressCoefficient::Eval(ElementTransformation &T,
380 const IntegrationPoint &ip)
381{
382 MFEM_ASSERT(u != NULL, "displacement field is not set");
383
384 real_t L = lambda.Eval(T, ip);
385 real_t M = mu.Eval(T, ip);
386 u->GetVectorGradient(T, grad);
387 if (si == sj)
388 {
389 real_t div_u = grad.Trace();
390 return L*div_u + 2*M*grad(si,si);
391 }
392 else
393 {
394 return M*(grad(si,sj) + grad(sj,si));
395 }
396}
397
398
399VisMan::VisMan(const char *vishost, const int visport)
400 : iostream(0),
401 host(vishost), port(visport), sid(0)
402{
403 win_x = 0;
404 win_y = 0;
405 win_w = 400; // window width
406 win_h = 350; // window height
407 win_stride_x = win_w;
408 win_stride_y = win_h + 20;
409 win_nx = 4; // number of windows in a row
410}
411
412void VisMan::NewWindow()
413{
414 sock.Append(new socketstream(host, port));
415 sid = sock.Size()-1;
416 iostream::rdbuf(sock[sid]->rdbuf());
417}
418
419void VisMan::CloseConnection()
420{
421 if (sid < sock.Size())
422 {
423 delete sock[sid];
424 sock[sid] = NULL;
425 iostream::rdbuf(0);
426 }
427}
428
429void VisMan::PositionWindow()
430{
431 *this << "window_geometry "
432 << win_x + win_stride_x*(sid%win_nx) << ' '
433 << win_y + win_stride_y*(sid/win_nx) << ' '
434 << win_w << ' ' << win_h << endl;
435}
436
437VisMan::~VisMan()
438{
439 for (int i = sock.Size()-1; i >= 0; i--)
440 {
441 delete sock[i];
442 }
443}
444
445ostream &operator<<(ostream &v, void (*f)(VisMan&))
446{
447 VisMan *vp = dynamic_cast<VisMan*>(&v);
448 if (vp) { (*f)(*vp); }
449 return v;
450}
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
@ GaussLobatto
Closed type.
Definition fe_base.hpp:32
Conjugate gradient method.
Definition solvers.hpp:513
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
real_t Trace() const
Trace of a square matrix.
Definition densemat.cpp:511
GMRES method.
Definition solvers.hpp:547
void SetKDim(int dim)
Set the number of iteration to perform between restarts, default is 50.
Definition solvers.hpp:559
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
void GetVectorGradient(ElementTransformation &tr, DenseMatrix &grad) const
Compute the vector gradient with respect to the physical element variable.
The BoomerAMG solver in hypre.
Definition hypre.hpp:1691
void SetSystemsOptions(int dim, bool order_bynodes=false)
Definition hypre.cpp:5111
void SetElasticityOptions(ParFiniteElementSpace *fespace, bool interp_refine=true)
Definition hypre.cpp:5238
Wrapper for hypre's ParCSR matrix class.
Definition hypre.hpp:388
MPI_Comm GetComm() const
MPI communicator.
Definition hypre.hpp:578
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition hypre.hpp:74
Class for integration point with weight.
Definition intrules.hpp:35
Abstract base class for iterative solver.
Definition solvers.hpp:67
virtual void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
Definition solvers.cpp:179
void SetRelTol(real_t rtol)
Definition solvers.hpp:209
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition solvers.cpp:173
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition solvers.cpp:71
void SetMaxIter(int max_it)
Definition solvers.hpp:211
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition fe_coll.hpp:330
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
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition mesh.hpp:290
void Clear()
Clear the contents of the Mesh.
Definition mesh.hpp:730
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1226
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
void GetNodes(Vector &node_coord) const
Definition mesh.cpp:8973
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 UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition mesh.hpp:280
static bool Root()
Return true if the rank in MPI_COMM_WORLD is zero.
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
static void Init(int &argc, char **&argv, int required=default_thread_required, int *provided=nullptr)
Singleton creation with Mpi::Init(argc, argv).
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
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
bool Good() const
Return true if the command line options were parsed successfully.
A piecewise constant coefficient with the constants keyed off the element attribute numbers.
Class for parallel bilinear form.
Abstract parallel finite element space.
Definition pfespace.hpp:29
HYPRE_BigInt GlobalTrueVSize() const
Definition pfespace.hpp:285
Class for parallel grid function.
Definition pgridfunc.hpp:33
void ProjectCoefficient(Coefficient &coeff) override
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Class for parallel linear form.
Class for parallel meshes.
Definition pmesh.hpp:34
int GetMyRank() const
Definition pmesh.hpp:404
int GetNRanks() const
Definition pmesh.hpp:403
void SetNodalFESpace(FiniteElementSpace *nfes) override
Definition pmesh.cpp:2028
A general vector function coefficient.
Vector data type.
Definition vector.hpp:80
void Neg()
(*this) = -(*this)
Definition vector.cpp:300
const real_t alpha
Definition ex15.cpp:369
void close_connection(VisMan &v)
Definition ex17p.cpp:98
void InitDisplacement(const Vector &x, Vector &u)
Definition ex17p.cpp:372
void position_window(VisMan &v)
Definition ex17p.cpp:97
void new_window(VisMan &v)
Definition ex17p.cpp:96
real_t kappa
Definition ex24.cpp:54
int dim
Definition ex24.cpp:53
real_t mu
Definition ex25.cpp:140
int main()
HYPRE_Int HYPRE_BigInt
real_t b
Definition lissajous.cpp:42
real_t a
Definition lissajous.cpp:41
std::ostream & operator<<(std::ostream &os, SparseMatrix const &mat)
const int visport
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
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[]