MFEM  v4.6.0
Finite element discretization library
ex2.cpp
Go to the documentation of this file.
1 // MFEM Example 2
2 // PUMI Modification
3 //
4 // Compile with: make ex2
5 //
6 // Sample runs:
7 // ex2 -m ../../data/pumi/serial/pillbox.smb -p ../../data/pumi/geom/pillbox.dmg
8 // -bf ../../data/pumi/serial/boundary.mesh
9 //
10 // Note: Example models + meshes for the PUMI examples can be downloaded
11 // from github.com/mfem/data/pumi. After downloading we recommend
12 // creating a symbolic link to the above directory in ../../data.
13 //
14 // Description: This example code solves a simple linear elasticity problem
15 // describing a multi-material cantilever beam.
16 //
17 // Specifically, we approximate the weak form of -div(sigma(u))=0
18 // where sigma(u)=lambda*div(u)*I+mu*(grad*u+u*grad) is the stress
19 // tensor corresponding to displacement field u, and lambda and mu
20 // are the material Lame constants. The boundary conditions are
21 // u=0 on the fixed part of the boundary with attribute 1, and
22 // sigma(u).n=f on the remainder with f being a constant pull down
23 // vector on boundary elements with attribute 2, and zero
24 // otherwise. The geometry of the domain is assumed to be as
25 // follows:
26 // boundary
27 // attribute 2
28 // (push down)
29 // ||
30 // \/
31 // +----------+
32 // | |
33 // | |
34 // +---------| material |----------+
35 // boundary --->| material| 2 | material |<--- boundary
36 // attribute 1 | 1 | | 3 | attribute 1
37 // (fixed) +---------+----------+----------+ (fixed)
38 //
39 // The example demonstrates the use of high-order and NURBS vector
40 // finite element spaces with the linear elasticity bilinear form,
41 // meshes with curved elements, and the definition of piece-wise
42 // constant and vector coefficient objects. Static condensation is
43 // also illustrated.
44 //
45 // We recommend viewing Example 1 before viewing this example.
46 //
47 // NOTE: Model/Mesh files for this example are in the (large) data file
48 // repository of MFEM here https://github.com/mfem/data under the
49 // folder named "pumi", which consists of the following sub-folders:
50 // a) geom --> model files
51 // b) parallel --> parallel pumi mesh files
52 // c) serial --> serial pumi mesh files
53 
54 
55 #include "mfem.hpp"
56 #include <fstream>
57 #include <iostream>
58 
59 #include "../../general/text.hpp"
60 
61 #ifdef MFEM_USE_SIMMETRIX
62 #include <SimUtil.h>
63 #include <gmi_sim.h>
64 #endif
65 #include <apfMDS.h>
66 #include <gmi_null.h>
67 #include <PCU.h>
68 #include <apfConvert.h>
69 #include <gmi_mesh.h>
70 #include <crv.h>
71 
72 #ifndef MFEM_USE_PUMI
73 #error This example requires that MFEM is built with MFEM_USE_PUMI=YES
74 #endif
75 
76 using namespace std;
77 using namespace mfem;
78 
79 int main(int argc, char *argv[])
80 {
81  // 1. Initialize MPI (required by PUMI) and HYPRE.
82  Mpi::Init(argc, argv);
83  int num_proc = Mpi::WorldSize();
84  int myId = Mpi::WorldRank();
85  Hypre::Init();
86 
87  // 2. Parse command-line options.
88  const char *mesh_file = "../../data/pumi/serial/pillbox.smb";
89  const char *boundary_file = "../../data/pumi/serial/boundary.mesh";
90 #ifdef MFEM_USE_SIMMETRIX
91  const char *model_file = "../../data/pumi/geom/pillbox.smd";
92 #else
93  const char *model_file = "../../data/pumi/geom/pillbox.dmg";
94 #endif
95  int order = 1;
96  bool static_cond = false;
97  bool visualization = 1;
98  int geom_order = 1;
99 
100  OptionsParser args(argc, argv);
101  args.AddOption(&mesh_file, "-m", "--mesh",
102  "Mesh file to use.");
103  args.AddOption(&order, "-o", "--order",
104  "Finite element order (polynomial degree).");
105  args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
106  "--no-static-condensation", "Enable static condensation.");
107  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
108  "--no-visualization",
109  "Enable or disable GLVis visualization.");
110  args.AddOption(&model_file, "-p", "--parasolid",
111  "Parasolid model to use.");
112  args.AddOption(&geom_order, "-go", "--geometry_order",
113  "Geometric order of the model");
114  args.AddOption(&boundary_file, "-bf", "--txt",
115  "txt file containing boundary tags");
116  args.Parse();
117  if (!args.Good())
118  {
119  args.PrintUsage(cout);
120  return 1;
121  }
122  args.PrintOptions(cout);
123 
124  // 3. Read the SCOREC Mesh.
125  PCU_Comm_Init();
126 #ifdef MFEM_USE_SIMMETRIX
127  Sim_readLicenseFile(0);
128  gmi_sim_start();
129  gmi_register_sim();
130 #endif
131  gmi_register_mesh();
132 
133  apf::Mesh2* pumi_mesh;
134  pumi_mesh = apf::loadMdsMesh(model_file, mesh_file);
135 
136  // 4. Increase the geometry order if necessary.
137  if (geom_order > 1)
138  {
139  crv::BezierCurver bc(pumi_mesh, geom_order, 0);
140  bc.run();
141  }
142  pumi_mesh->verify();
143 
144  // Read boundary
145  string bdr_tags;
146  named_ifgzstream input_bdr(boundary_file);
147  input_bdr >> ws;
148  getline(input_bdr, bdr_tags);
149  filter_dos(bdr_tags);
150  cout << " the boundary tag is : " << bdr_tags << endl;
151  Array<int> Dirichlet;
152  int numOfent;
153  if (bdr_tags == "Dirichlet")
154  {
155  input_bdr >> numOfent;
156  cout << " num of Dirichlet bdr conditions : " << numOfent << endl;
157  Dirichlet.SetSize(numOfent);
158  for (int kk = 0; kk < numOfent; kk++)
159  {
160  input_bdr >> Dirichlet[kk];
161  }
162  }
163  Dirichlet.Print();
164 
165  Array<int> load_bdr;
166  skip_comment_lines(input_bdr, '#');
167  input_bdr >> bdr_tags;
168  filter_dos(bdr_tags);
169  cout << " the boundary tag is : " << bdr_tags << endl;
170  if (bdr_tags == "Load")
171  {
172  input_bdr >> numOfent;
173  load_bdr.SetSize(numOfent);
174  cout << " num of load bdr conditions : " << numOfent << endl;
175  for (int kk = 0; kk < numOfent; kk++)
176  {
177  input_bdr >> load_bdr[kk];
178  }
179  }
180  load_bdr.Print();
181 
182  // 5. Create the MFEM mesh object from the PUMI mesh. We can handle triangular
183  // and tetrahedral meshes. Other inputs are the same as MFEM default
184  // constructor.
185  Mesh *mesh = new PumiMesh(pumi_mesh, 1, 1);
186  int dim = mesh->Dimension();
187 
188  // Boundary conditions hack.
189  apf::MeshIterator* itr = pumi_mesh->begin(dim-1);
190  apf::MeshEntity* ent ;
191  int bdr_cnt = 0;
192  while ((ent = pumi_mesh->iterate(itr)))
193  {
194  apf::ModelEntity *me = pumi_mesh->toModel(ent);
195  if (pumi_mesh->getModelType(me) == (dim-1))
196  {
197  // Everywhere 3 as initial
198  (mesh->GetBdrElement(bdr_cnt))->SetAttribute(3);
199  int tag = pumi_mesh->getModelTag(me);
200  if (Dirichlet.Find(tag) != -1)
201  {
202  // Dirichlet attr -> 1
203  (mesh->GetBdrElement(bdr_cnt))->SetAttribute(1);
204  }
205  else if (load_bdr.Find(tag) != -1)
206  {
207  // Load attr -> 2
208  (mesh->GetBdrElement(bdr_cnt))->SetAttribute(2);
209  }
210  bdr_cnt++;
211  }
212  }
213  pumi_mesh->end(itr);
214 
215  // Assign attributes for elements.
216  double ppt[3];
217  Vector cent(ppt, dim);
218  for (int el = 0; el < mesh->GetNE(); el++)
219  {
220  (mesh->GetElementTransformation(el))->
221  Transform(Geometries.GetCenter(mesh->GetElementBaseGeometry(el)),cent);
222  if (cent(0) <= -0.05)
223  {
224  mesh->SetAttribute(el, 1);
225  }
226  else if (cent(0) >= 0.05)
227  {
228  mesh->SetAttribute(el, 2);
229  }
230  else
231  {
232  mesh->SetAttribute(el, 3);
233  }
234  }
235  mesh->SetAttributes();
236  if (mesh->attributes.Max() < 2 || mesh->bdr_attributes.Max() < 2)
237  {
238  cerr << "\nInput mesh should have at least two materials and "
239  << "two boundary attributes! (See schematic in ex2.cpp)\n"
240  << endl;
241  return 3;
242  }
243 
244  // 6. Refine the mesh to increase the resolution. In this example we do
245  // 'ref_levels' of uniform refinement. We choose 'ref_levels' to be the
246  // largest number that gives a final mesh with no more than 5,000
247  // elements.
248  {
249  int ref_levels =
250  (int)floor(log(5000./mesh->GetNE())/log(2.)/dim);
251  for (int l = 0; l < ref_levels; l++)
252  {
253  mesh->UniformRefinement();
254  }
255  }
256 
257  // 7. Define a finite element space on the mesh. Here we use vector finite
258  // elements, i.e. dim copies of a scalar finite element space. The vector
259  // dimension is specified by the last argument of the FiniteElementSpace
260  // constructor. For NURBS meshes, we use the (degree elevated) NURBS space
261  // associated with the mesh nodes.
263  FiniteElementSpace *fespace;
264  if (mesh->NURBSext)
265  {
266  fec = NULL;
267  fespace = mesh->GetNodes()->FESpace();
268  }
269  else
270  {
271  fec = new H1_FECollection(order, dim);
272  fespace = new FiniteElementSpace(mesh, fec, dim);
273  }
274  cout << "Number of finite element unknowns: " << fespace->GetTrueVSize()
275  << endl << "Assembling: " << flush;
276 
277  // 8. Determine the list of true (i.e. conforming) essential boundary dofs.
278  // In this example, the boundary conditions are defined by marking only
279  // boundary attribute 1 from the mesh as essential and converting it to a
280  // list of true dofs.
281  Array<int> ess_tdof_list, ess_bdr(mesh->bdr_attributes.Max());
282  ess_bdr = 0;
283  ess_bdr[0] = 1;
284  fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
285 
286  // 9. Set up the linear form b(.) which corresponds to the right-hand side of
287  // the FEM linear system. In this case, b_i equals the boundary integral
288  // of f*phi_i where f represents a "pull down" force on the Neumann part
289  // of the boundary and phi_i are the basis functions in the finite element
290  // fespace. The force is defined by the VectorArrayCoefficient object f,
291  // which is a vector of Coefficient objects. The fact that f is non-zero
292  // on boundary attribute 2 is indicated by the use of piece-wise constants
293  // coefficient for its last component.
295  for (int i = 0; i < dim-1; i++)
296  {
297  f.Set(i, new ConstantCoefficient(0.0));
298  }
299  {
300  Vector pull_force(mesh->bdr_attributes.Max());
301  pull_force = 0.0;
302  pull_force(1) = -3.0e-2;
303  f.Set(dim-1, new PWConstCoefficient(pull_force));
304  f.Set(dim-2, new PWConstCoefficient(pull_force));
305  }
306 
307  LinearForm *b = new LinearForm(fespace);
308  b->AddBoundaryIntegrator(new VectorBoundaryLFIntegrator(f));
309  cout << "r.h.s. ... " << flush;
310  b->Assemble();
311 
312  // 10. Define the solution vector x as a finite element grid function
313  // corresponding to fespace. Initialize x with initial guess of zero,
314  // which satisfies the boundary conditions.
315  GridFunction x(fespace);
316  x = 0.0;
317 
318  // 11. Set up the bilinear form a(.,.) on the finite element space
319  // corresponding to the linear elasticity integrator with piece-wise
320  // constants coefficient lambda and mu.
321  Vector lambda(mesh->attributes.Max());
322  lambda = 1.0;
323  lambda(0) = lambda(1)*10;
324  lambda(1) = lambda(1)*100;
325  PWConstCoefficient lambda_func(lambda);
326  Vector mu(mesh->attributes.Max());
327  mu = 1.0;
328  mu(0) = mu(1)*10;
329  mu(1) = mu(1)*100;
330  PWConstCoefficient mu_func(mu);
331 
332  BilinearForm *a = new BilinearForm(fespace);
333  a->AddDomainIntegrator(new ElasticityIntegrator(lambda_func,mu_func));
334 
335  // 12. Assemble the bilinear form and the corresponding linear system,
336  // applying any necessary transformations such as: eliminating boundary
337  // conditions, applying conforming constraints for non-conforming AMR,
338  // static condensation, etc.
339  cout << "matrix ... " << flush;
340  if (static_cond) { a->EnableStaticCondensation(); }
341  a->Assemble();
342 
343  SparseMatrix A;
344  Vector B, X;
345  a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
346  cout << "done." << endl;
347 
348  cout << "Size of linear system: " << A.Height() << endl;
349 
350 #ifndef MFEM_USE_SUITESPARSE
351  // 13. Define a simple symmetric Gauss-Seidel preconditioner and use it to
352  // solve the system Ax=b with PCG.
353  GSSmoother M(A);
354  PCG(A, M, B, X, 1, 500, 1e-8, 0.0);
355 #else
356  // 13. If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
357  UMFPackSolver umf_solver;
358  umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
359  umf_solver.SetOperator(A);
360  umf_solver.Mult(B, X);
361 #endif
362 
363  // 14. Recover the solution as a finite element grid function.
364  a->RecoverFEMSolution(X, *b, x);
365 
366  // 15. For non-NURBS meshes, make the mesh curved based on the finite element
367  // space. This means that we define the mesh elements through a fespace
368  // based transformation of the reference element. This allows us to save
369  // the displaced mesh as a curved mesh when using high-order finite
370  // element displacement field. We assume that the initial mesh (read from
371  // the file) is not higher order curved mesh compared to the chosen FE
372  // space.
373  if (!mesh->NURBSext)
374  {
375  mesh->SetNodalFESpace(fespace);
376  }
377 
378  // 16. Save the displaced mesh and the inverted solution (which gives the
379  // backward displacements to the original grid). This output can be
380  // viewed later using GLVis: "glvis -m displaced.mesh -g sol.gf".
381  {
382  GridFunction *nodes = mesh->GetNodes();
383  *nodes += x;
384  x *= -1;
385  ofstream mesh_ofs("displaced.mesh");
386  mesh_ofs.precision(8);
387  mesh->Print(mesh_ofs);
388  ofstream sol_ofs("sol.gf");
389  sol_ofs.precision(8);
390  x.Save(sol_ofs);
391  }
392 
393  // 17. Send the above data by socket to a GLVis server. Use the "n" and "b"
394  // keys in GLVis to visualize the displacements.
395  if (visualization)
396  {
397  char vishost[] = "localhost";
398  int visport = 19916;
399  socketstream sol_sock(vishost, visport);
400  sol_sock.precision(8);
401  sol_sock << "solution\n" << *mesh << x << flush;
402  }
403 
404  // 18. Free the used memory.
405  delete a;
406  delete b;
407  if (fec)
408  {
409  delete fespace;
410  delete fec;
411  }
412  delete mesh;
413 
414  pumi_mesh->destroyNative();
415  apf::destroyMesh(pumi_mesh);
416  PCU_Comm_Free();
417 #ifdef MFEM_USE_SIMMETRIX
418  gmi_sim_stop();
419  Sim_unregisterAllKeys();
420 #endif
421 
422  return 0;
423 }
int visport
Vector coefficient defined by an array of scalar coefficients. Coefficients that are not set will eva...
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:331
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:1242
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:68
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
Base class for PUMI meshes.
Definition: pumi.hpp:44
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:159
STL namespace.
const IntegrationPoint & GetCenter(int GeomType)
Return the center of the given Geometry::Type, GeomType.
Definition: geom.hpp:71
Geometry Geometries
Definition: fe.cpp:49
void skip_comment_lines(std::istream &is, const char comment_char)
Check if the stream starts with comment_char. If so skip it.
Definition: text.hpp:31
Data type for Gauss-Seidel smoother of sparse matrix.
Direct sparse solver using UMFPACK.
Definition: solvers.hpp:1070
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
char vishost[]
Data type sparse matrix.
Definition: sparsemat.hpp:50
virtual void SetAttributes()
Determine the sets of unique attribute values in domain and boundary elements.
Definition: mesh.cpp:1572
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:10232
int Find(const T &el) const
Return the first index where &#39;el&#39; is found; return -1 if not found.
Definition: array.hpp:818
void Print(std::ostream &out=mfem::out, int width=4) const
Prints array to stream with width elements per row.
Definition: array.cpp:23
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
virtual void SetNodalFESpace(FiniteElementSpace *nfes)
Definition: mesh.cpp:5577
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:712
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:275
double Control[UMFPACK_CONTROL]
Definition: solvers.hpp:1081
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:219
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
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
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:687
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:3194
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1086
double a
Definition: lissajous.cpp:41
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:277
A piecewise constant coefficient with the constants keyed off the element attribute numbers...
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
int main(int argc, char *argv[])
Definition: ex2.cpp:47
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: mesh.cpp:355
double mu
Definition: ex25.cpp:140
int dim
Definition: ex24.cpp:53
void filter_dos(std::string &line)
Check for, and remove, a trailing &#39;\r&#39; from and std::string.
Definition: text.hpp:45
void SetAttribute(int i, int attr)
Set the attribute of element i.
Definition: mesh.hpp:1193
Vector data type.
Definition: vector.hpp:58
virtual void Print(std::ostream &os=mfem::out) const
Definition: mesh.hpp:2011
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:259
Vector with associated FE space and LinearFormIntegrators.
Definition: linearform.hpp:24
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:8302
const Element * GetBdrElement(int i) const
Return pointer to the i&#39;th boundary element object.
Definition: mesh.hpp:1158
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:3696
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition: mesh.hpp:273
double f(const Vector &p)
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a SparseMatrix.
Definition: solvers.cpp:3099