MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
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
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
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
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
76using namespace std;
77using namespace mfem;
78
79int 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();
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);
102 "Mesh file to use.");
104 "Finite element order (polynomial degree).");
106 "--no-static-condensation", "Enable static condensation.");
108 "--no-visualization",
109 "Enable or disable GLVis visualization.");
111 "Parasolid model to use.");
113 "Geometric order of the model");
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
128 gmi_sim_start();
129 gmi_register_sim();
130#endif
131 gmi_register_mesh();
132
133 apf::Mesh2* pumi_mesh;
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
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
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;
171 {
172 input_bdr >> numOfent;
174 cout << " num of load bdr conditions : " << numOfent << endl;
175 for (int kk = 0; kk < numOfent; kk++)
176 {
178 }
179 }
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);
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);
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}
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition array.hpp:697
int Find(const T &el) const
Return the first index where 'el' is found; return -1 if not found.
Definition array.hpp:828
void Print(std::ostream &out=mfem::out, int width=4) const
Prints array to stream with width elements per row.
Definition array.cpp:23
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.
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition fe_coll.hpp:27
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
Data type for Gauss-Seidel smoother of sparse matrix.
const IntegrationPoint & GetCenter(int GeomType) const
Return the center of the given Geometry::Type, GeomType.
Definition geom.hpp:71
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.
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition hypre.hpp:74
Vector with associated FE space and LinearFormIntegrators.
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 SetAttribute(int i, int attr)
Set the attribute of element i.
Definition mesh.hpp:1336
virtual void Print(std::ostream &os=mfem::out, const std::string &comments="") const
Definition mesh.hpp:2288
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
const Element * GetBdrElement(int i) const
Return pointer to the i'th boundary element object.
Definition mesh.hpp:1298
void GetElementTransformation(int i, IsoparametricTransformation *ElTr) const
Builds the transformation defining the i-th element in ElTr. ElTr must be allocated in advance and wi...
Definition mesh.cpp:357
void GetNodes(Vector &node_coord) const
Definition mesh.cpp:8973
virtual void SetNodalFESpace(FiniteElementSpace *nfes)
Definition mesh.cpp:6153
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
Geometry::Type GetElementBaseGeometry(int i) const
Definition mesh.hpp:1385
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition mesh.hpp:280
virtual void SetAttributes()
Determine the sets of unique attribute values in domain and boundary elements.
Definition mesh.cpp:1604
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
static int WorldSize()
Return the size of 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).
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition operator.hpp:66
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.
Base class for PUMI meshes.
Definition pumi.hpp:45
Data type sparse matrix.
Definition sparsemat.hpp:51
Direct sparse solver using UMFPACK.
Definition solvers.hpp:1096
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a SparseMatrix.
Definition solvers.cpp:3102
real_t Control[UMFPACK_CONTROL]
Definition solvers.hpp:1106
virtual void Mult(const Vector &b, Vector &x) const
Direct solution of the linear system using UMFPACK.
Definition solvers.cpp:3197
Vector coefficient defined by an array of scalar coefficients. Coefficients that are not set will eva...
Vector data type.
Definition vector.hpp:80
Vector & Set(const real_t a, const Vector &x)
(*this) = a * x
Definition vector.cpp:262
int dim
Definition ex24.cpp:53
real_t mu
Definition ex25.cpp:140
int main()
real_t b
Definition lissajous.cpp:42
real_t a
Definition lissajous.cpp:41
const int visport
Geometry Geometries
Definition fe.cpp:49
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
void filter_dos(std::string &line)
Check for, and remove, a trailing '\r' from and std::string.
Definition text.hpp:45
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30
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
const char vishost[]