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