MFEM  v4.5.2
Finite element discretization library
display-basis.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 //
12 // ----------------------------------------------------------------
13 // Display Basis Miniapp: Visualize finite element basis functions
14 // ----------------------------------------------------------------
15 //
16 // This miniapp visualizes various types of finite element basis functions on a
17 // single mesh element in 1D, 2D and 3D. The order and the type of finite
18 // element space can be changed, and the mesh element is either the reference
19 // one, or a simple transformation of it. Dynamic creation and interaction with
20 // multiple GLVis windows is demonstrated.
21 //
22 // Compile with: make display-basis
23 //
24 // Sample runs: display-basis
25 // display_basis -e 2 -b 3 -o 3
26 // display-basis -e 5 -b 1 -o 1
27 // display-basis -e 3 -b 7 -o 3
28 // display-basis -e 3 -b 7 -o 5 -only 16
29 
30 #include "mfem.hpp"
31 #include "../common/mfem-common.hpp"
32 #include <vector>
33 #include <iostream>
34 
35 using namespace std;
36 using namespace mfem;
37 using namespace mfem::common;
38 
39 // Data structure used to collect visualization window layout parameters
40 struct VisWinLayout
41 {
42  int nx;
43  int ny;
44  int w;
45  int h;
46 };
47 
48 // Data structure used to define simple coordinate transformations
49 struct DeformationData
50 {
51  double uniformScale;
52 
53  int squeezeAxis;
54  double squeezeFactor;
55 
56  int shearAxis;
57  Vector shearVec;
58 };
59 
60 /** The Deformation class implements three simple coordinate transformations:
61  Uniform Scaling:
62  u = a v for a scalar constant 'a'
63 
64  Compression or Squeeze (along a coordinate axis):
65  / 1/b 0 \ / 1/b 0 0 \ for a scalar constant b
66  u = \ 0 b / v or u = | 0 c 0 | v, and c = sqrt(b)
67  \ 0 0 c / the axis can also be chosen
68 
69  Shear:
70  u = v + v_i * s where 's' is the shear vector
71  and 'i' is the shear axis
72 */
73 class Deformation : public VectorCoefficient
74 {
75 public:
76 
77  enum DefType {INVALID, UNIFORM, SQUEEZE, SHEAR};
78 
79  Deformation(int dim, DefType dType, const DeformationData & data)
80  : VectorCoefficient(dim), dim_(dim), dType_(dType), data_(data) {}
81 
82  void Eval(Vector &v, ElementTransformation &T, const IntegrationPoint &ip);
83  using VectorCoefficient::Eval;
84 private:
85  void Def1D(const Vector & u, Vector & v);
86  void Def2D(const Vector & u, Vector & v);
87  void Def3D(const Vector & u, Vector & v);
88 
89  int dim_;
90  DefType dType_;
91  const DeformationData & data_;
92 };
93 
94 string elemTypeStr(const Element::Type & eType);
95 inline bool elemIs1D(const Element::Type & eType);
96 inline bool elemIs2D(const Element::Type & eType);
97 inline bool elemIs3D(const Element::Type & eType);
98 
99 string basisTypeStr(char bType);
100 inline bool basisIs1D(char bType);
101 inline bool basisIs2D(char bType);
102 inline bool basisIs3D(char bType);
103 
104 string mapTypeStr(int mType);
105 
106 int update_basis(vector<socketstream*> & sock, const VisWinLayout & vwl,
107  Element::Type e, char bType, int bOrder, int mType,
108  Deformation::DefType dType, const DeformationData & defData,
109  bool visualization, int &onlySome);
110 
111 int main(int argc, char *argv[])
112 {
113  // Parse command-line options.
114  Element::Type eType = Element::TRIANGLE;
115  char bType = 'h';
116  int bOrder = 2;
117  int mType = 0;
118 
119  int eInt = -1;
120  int bInt = -1;
121 
122  VisWinLayout vwl;
123  vwl.nx = 5;
124  vwl.ny = 3;
125  vwl.w = 250;
126  vwl.h = 250;
127 
128  Deformation::DefType dType = Deformation::INVALID;
129  DeformationData defData;
130 
131  bool visualization = true;
132  int onlySome = -1;
133 
134  vector<socketstream*> sock;
135 
136  OptionsParser args(argc, argv);
137  args.AddOption(&eInt, "-e", "--elem-type",
138  "Element Type: (1-Segment, 2-Triangle, 3-Quadrilateral, "
139  "4-Tetrahedron, 5-Hexahedron, 6-Wedge)");
140  args.AddOption(&bInt, "-b", "--basis-type",
141  "Basis Function Type (0-H1, 1-Nedelec, 2-Raviart-Thomas, "
142  "3-L2, 4-Fixed Order Cont.,\n\t5-Gaussian Discontinuous (2D),"
143  " 6-Crouzeix-Raviart, 7-Serendipity)");
144  args.AddOption(&bOrder, "-o", "--order", "Basis function order");
145  args.AddOption(&vwl.nx, "-nx", "--num-win-x",
146  "Number of Viz windows in X");
147  args.AddOption(&vwl.ny, "-ny", "--num-win-y",
148  "Number of Viz windows in y");
149  args.AddOption(&vwl.w, "-w", "--width",
150  "Width of Viz windows");
151  args.AddOption(&vwl.h, "-h", "--height",
152  "Height of Viz windows");
153  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
154  "--no-visualization",
155  "Enable or disable GLVis visualization.");
156  args.AddOption(&onlySome, "-only", "--onlySome",
157  "Only view 10 dofs, starting with the specified one.");
158  args.Parse();
159  if (!args.Good())
160  {
161  args.PrintUsage(cout);
162  return 1;
163  }
164  {
165  args.PrintOptions(cout);
166  }
167  if ( eInt > 0 && eInt < 7 )
168  {
169  eType = (Element::Type)eInt;
170  }
171  switch (bInt)
172  {
173  case 0:
174  bType = 'h';
175  break;
176  case 1:
177  bType = 'n';
178  break;
179  case 2:
180  bType = 'r';
181  break;
182  case 3:
183  bType = 'l';
184  break;
185  case 4:
186  bType = 'f';
187  break;
188  case 5:
189  bType = 'g';
190  break;
191  case 6:
192  bType = 'c';
193  break;
194  case 7:
195  bType = 's';
196  break;
197  default:
198  bType = 'h';
199  }
200 
201  // Collect user input
202  bool print_char = true;
203  while (true)
204  {
205  if (print_char)
206  {
207  cout << endl;
208  cout << "Element Type: " << elemTypeStr(eType) << endl;
209  cout << "Basis Type: " << basisTypeStr(bType) << endl;;
210  cout << "Basis function order: " << bOrder << endl;
211  cout << "Map Type: " << mapTypeStr(mType) << endl;
212  }
213  if ( update_basis(sock, vwl, eType, bType, bOrder, mType,
214  dType, defData, visualization, onlySome) )
215  {
216  cerr << "Invalid combination of basis info (try again)" << endl;
217  }
218 
219  if (!visualization) { break; }
220 
221  print_char = false;
222  cout << endl;
223  cout << "What would you like to do?\n"
224  "q) Quit\n"
225  "c) Close Windows and Quit\n"
226  "e) Change Element Type\n"
227  "b) Change Basis Type\n";
228  if ( bType == 'h' || bType == 'p' || bType == 'n' || bType == 'r' ||
229  bType == 'l' || bType == 'f' || bType == 'g' || bType == 's')
230  {
231  cout << "o) Change Basis Order\n";
232  }
233  // The following is disabled pending updates to GLVis
234  if ( bType == 'l' && false )
235  {
236  cout << "m) Change Map Type\n";
237  }
238  cout << "t) Transform Element\n";
239  cout << "--> " << flush;
240  char mk;
241  cin >> mk;
242 
243  if (mk == 'q')
244  {
245  break;
246  }
247  if (mk == 'c')
248  {
249  for (unsigned int i=0; i<sock.size(); i++)
250  {
251  *sock[i] << "keys q";
252  }
253  break;
254  }
255  if (mk == 'e')
256  {
257  eInt = 0;
258  cout << "valid element types:\n";
259  if ( basisIs1D(bType) )
260  {
261  cout <<
262  "1) Segment\n";
263  }
264  if ( basisIs2D(bType) )
265  {
266  cout <<
267  "2) Triangle\n"
268  "3) Quadrilateral\n";
269  }
270  if ( basisIs3D(bType) )
271  {
272  cout <<
273  "4) Tetrahedron\n"
274  "5) Hexahedron\n"
275  "6) Wedge\n";
276  }
277  cout << "enter new element type --> " << flush;
278  cin >> eInt;
279  if ( eInt <= 0 || eInt > 6 )
280  {
281  cout << "invalid element type \"" << eInt << "\"" << endl << flush;
282  }
283  else if ( (elemIs1D((Element::Type)eInt) && basisIs1D(bType)) ||
284  (elemIs2D((Element::Type)eInt) && basisIs2D(bType)) ||
285  (elemIs3D((Element::Type)eInt) && basisIs3D(bType)) )
286  {
287  if ( (elemIs1D((Element::Type)eInt) && !elemIs1D(eType)) ||
288  (elemIs2D((Element::Type)eInt) && !elemIs2D(eType)) ||
289  (elemIs3D((Element::Type)eInt) && !elemIs3D(eType)) )
290  {
291  dType = Deformation::INVALID;
292  }
293  eType = (Element::Type)eInt;
294 
295  print_char = true;
296  }
297  else
298  {
299  cout << "invalid element type \"" << eInt <<
300  "\" for basis type \"" << basisTypeStr(bType) << "\"." << endl;
301  }
302  }
303  if (mk == 'b')
304  {
305  char bChar = 0;
306  cout << "valid basis types:\n";
307  cout << "h) H1 Finite Element\n";
308  cout << "p) H1 Positive Finite Element\n";
309  if ( elemIs2D(eType) || elemIs3D(eType) )
310  {
311  cout << "s) H1 Serendipity Finite Element\n";
312  cout << "n) Nedelec Finite Element\n";
313  cout << "r) Raviart-Thomas Finite Element\n";
314  }
315  cout << "l) L2 Finite Element\n";
316  if ( elemIs1D(eType) || elemIs2D(eType) )
317  {
318  cout << "c) Crouzeix-Raviart Finite Element\n";
319  }
320  cout << "f) Fixed Order Continuous Finite Element\n";
321  if ( elemIs2D(eType) )
322  {
323  cout << "g) Gauss Discontinuous Finite Element\n";
324  }
325  cout << "enter new basis type --> " << flush;
326  cin >> bChar;
327  if (bChar == 'h' || bChar == 'p' || bChar == 'l' || bChar == 'f' ||
328  bChar == 's' ||
329  ((bChar == 'n' || bChar == 'r') && (elemIs2D(eType) || elemIs3D(eType))) ||
330  (bChar == 'c' && (elemIs1D(eType) || elemIs2D(eType))) ||
331  (bChar == 'g' && elemIs2D(eType)))
332  {
333  bType = bChar;
334  if ( bType == 'h' )
335  {
336  mType = FiniteElement::VALUE;
337  }
338  else if ( bType == 'p' )
339  {
340  mType = FiniteElement::VALUE;
341  }
342  else if (bType == 's')
343  {
344  mType = FiniteElement::VALUE;
345  }
346  else if ( bType == 'n' )
347  {
348  mType = FiniteElement::H_CURL;
349  }
350  else if ( bType == 'r' )
351  {
352  mType = FiniteElement::H_DIV;
353  }
354  else if ( bType == 'l' )
355  {
356  if ( mType != FiniteElement::VALUE &&
357  mType != FiniteElement::INTEGRAL )
358  {
359  mType = FiniteElement::VALUE;
360  }
361  }
362  else if ( bType == 'c' )
363  {
364  bOrder = 1;
365  mType = FiniteElement::VALUE;
366  }
367  else if ( bType == 'f' )
368  {
369  if ( bOrder < 1 || bOrder > 3)
370  {
371  bOrder = 1;
372  }
373  mType = FiniteElement::VALUE;
374  }
375  else if ( bType == 'g' )
376  {
377  if ( bOrder < 1 || bOrder > 2)
378  {
379  bOrder = 1;
380  }
381  mType = FiniteElement::VALUE;
382  }
383  print_char = true;
384  }
385  else
386  {
387  cout << "invalid basis type \"" << bChar << "\"." << endl;
388  }
389  }
390  if (mk == 'm' && bType == 'l')
391  {
392  int mInt = 0;
393  cout << "valid map types:\n"
394  "0) VALUE\n"
395  "1) INTEGRAL\n";
396  cout << "enter new map type --> " << flush;
397  cin >> mInt;
398  if (mInt >=0 && mInt <= 1)
399  {
400  mType = mInt;
401  print_char = true;
402  }
403  else
404  {
405  cout << "invalid map type \"" << mInt << "\"." << endl;
406  }
407  }
408  if (mk == 'o')
409  {
410  int oInt = 1;
411  int oMin = ( bType == 'h' || bType == 'p' || bType == 'n' ||
412  bType == 'f' || bType == 'g' || bType == 's')?1:0;
413  int oMax = -1;
414  switch (bType)
415  {
416  case 'g':
417  oMax = 2;
418  break;
419  case 'f':
420  oMax = 3;
421  break;
422  default:
423  oMax = -1;
424  }
425  cout << "basis function order must be >= " << oMin;
426  if ( oMax >= 0 )
427  {
428  cout << " and <= " << oMax;
429  }
430  cout << endl;
431  cout << "enter new basis function order --> " << flush;
432  cin >> oInt;
433  if ( oInt >= oMin && oInt <= (oMax>=0)?oMax:oInt )
434  {
435  bOrder = oInt;
436  print_char = true;
437  }
438  else
439  {
440  cout << "invalid basis order \"" << oInt << "\"." << endl;
441  }
442  }
443  if (mk == 't')
444  {
445  cout << "transformation options:\n";
446  cout << "r) reset to reference element\n";
447  cout << "u) uniform scaling\n";
448  if ( elemIs2D(eType) || elemIs3D(eType) )
449  {
450  cout << "c) compression\n";
451  cout << "s) shear\n";
452  }
453  cout << "enter transformation type --> " << flush;
454  char tk;
455  cin >> tk;
456  if (tk == 'r')
457  {
458  dType = Deformation::INVALID;
459  }
460  else if (tk == 'u')
461  {
462  cout << "enter scaling constant --> " << flush;
463  cin >> defData.uniformScale;
464  if ( defData.uniformScale > 0.0 )
465  {
466  dType = Deformation::UNIFORM;
467  }
468  }
469  else if (tk == 'c' && !elemIs1D(eType))
470  {
471  int dim = elemIs2D(eType)?2:3;
472  cout << "enter compression factor --> " << flush;
473  cin >> defData.squeezeFactor;
474  cout << "enter compression axis (0-" << dim-1 << ") --> " << flush;
475  cin >> defData.squeezeAxis;
476 
477  if ( defData.squeezeFactor > 0.0 &&
478  (defData.squeezeAxis >= 0 && defData.squeezeAxis < dim))
479  {
480  dType = Deformation::SQUEEZE;
481  }
482  }
483  else if (tk == 's' && !elemIs1D(eType))
484  {
485  int dim = elemIs2D(eType)?2:3;
486  cout << "enter shear vector (components separated by spaces) --> "
487  << flush;
488  defData.shearVec.SetSize(dim);
489  for (int i=0; i<dim; i++)
490  {
491  cin >> defData.shearVec[i];
492  }
493  cout << "enter shear axis (0-" << dim-1 << ") --> " << flush;
494  cin >> defData.shearAxis;
495 
496  if ( defData.shearAxis >= 0 && defData.shearAxis < dim )
497  {
498  dType = Deformation::SHEAR;
499  }
500  }
501  }
502  }
503 
504  // Cleanup
505  for (unsigned int i=0; i<sock.size(); i++)
506  {
507  delete sock[i];
508  }
509 
510  // Exit
511  return 0;
512 }
513 
514 string elemTypeStr(const Element::Type & eType)
515 {
516  switch (eType)
517  {
518  case Element::POINT:
519  return "POINT";
520  case Element::SEGMENT:
521  return "SEGMENT";
522  case Element::TRIANGLE:
523  return "TRIANGLE";
524  case Element::QUADRILATERAL:
525  return "QUADRILATERAL";
526  case Element::TETRAHEDRON:
527  return "TETRAHEDRON";
528  case Element::HEXAHEDRON:
529  return "HEXAHEDRON";
530  case Element::WEDGE:
531  return "WEDGE";
532  default:
533  return "INVALID";
534  };
535 }
536 
537 bool
538 elemIs1D(const Element::Type & eType)
539 {
540  return eType == Element::SEGMENT;
541 }
542 
543 bool
544 elemIs2D(const Element::Type & eType)
545 {
546  return eType == Element::TRIANGLE || eType == Element::QUADRILATERAL;
547 }
548 
549 bool
550 elemIs3D(const Element::Type & eType)
551 {
552  return eType == Element::TETRAHEDRON || eType == Element::HEXAHEDRON ||
553  eType == Element::WEDGE;
554 }
555 
556 string
557 basisTypeStr(char bType)
558 {
559  switch (bType)
560  {
561  case 'h':
562  return "Continuous (H1)";
563  case 'p':
564  return "Continuous Positive (H1)";
565  case 's':
566  return "Continuous Serendipity (H1)";
567  case 'n':
568  return "Nedelec";
569  case 'r':
570  return "Raviart-Thomas";
571  case 'l':
572  return "Discontinuous (L2)";
573  case 'f':
574  return "Fixed Order Continuous";
575  case 'g':
576  return "Gaussian Discontinuous";
577  case 'c':
578  return "Crouzeix-Raviart";
579  default:
580  return "INVALID";
581  };
582 }
583 
584 bool
585 basisIs1D(char bType)
586 {
587  return bType == 'h' || bType == 'p' || bType == 'l' || bType == 'c' ||
588  bType == 'f';
589 }
590 
591 bool
592 basisIs2D(char bType)
593 {
594  return bType == 'h' || bType == 'p' || bType == 'n' || bType == 'r' ||
595  bType == 'l' || bType == 'c' || bType == 'f' || bType == 'g' ||
596  bType == 's';
597 }
598 
599 bool
600 basisIs3D(char bType)
601 {
602  return bType == 'h' || bType == 'p' || bType == 'n' || bType == 'r' ||
603  bType == 'f' || bType == 'l';
604 }
605 
606 string
607 mapTypeStr(int mType)
608 {
609  switch (mType)
610  {
611  case FiniteElement::VALUE:
612  return "VALUE";
613  case FiniteElement::H_CURL:
614  return "H_CURL";
615  case FiniteElement::H_DIV:
616  return "H_DIV";
617  case FiniteElement::INTEGRAL:
618  return "INTEGRAL";
619  default:
620  return "INVALID";
621  }
622 }
623 
624 void
626  const IntegrationPoint &ip)
627 {
628  Vector u(dim_);
629  T.Transform(ip, u);
630 
631  switch (dim_)
632  {
633  case 1:
634  Def1D(u, v);
635  break;
636  case 2:
637  Def2D(u, v);
638  break;
639  case 3:
640  Def3D(u, v);
641  break;
642  }
643 }
644 
645 void
646 Deformation::Def1D(const Vector & u, Vector & v)
647 {
648  v = u;
649  if ( dType_ == UNIFORM )
650  {
651  v *= data_.uniformScale;
652  }
653 }
654 
655 void
656 Deformation::Def2D(const Vector & u, Vector & v)
657 {
658  switch (dType_)
659  {
660  case UNIFORM:
661  v = u;
662  v *= data_.uniformScale;
663  break;
664  case SQUEEZE:
665  v = u;
666  v[ data_.squeezeAxis ] /= data_.squeezeFactor;
667  v[(data_.squeezeAxis+1)%2] *= data_.squeezeFactor;
668  break;
669  case SHEAR:
670  v = u;
671  v.Add(v[data_.shearAxis], data_.shearVec);
672  break;
673  default:
674  v = u;
675  }
676 }
677 
678 void
679 Deformation::Def3D(const Vector & u, Vector & v)
680 {
681  switch (dType_)
682  {
683  case UNIFORM:
684  v = u;
685  v *= data_.uniformScale;
686  break;
687  case SQUEEZE:
688  v = u;
689  v[ data_.squeezeAxis ] /= data_.squeezeFactor;
690  v[(data_.squeezeAxis+1)%2] *= sqrt(data_.squeezeFactor);
691  v[(data_.squeezeAxis+2)%2] *= sqrt(data_.squeezeFactor);
692  break;
693  case SHEAR:
694  v = u;
695  v.Add(v[data_.shearAxis], data_.shearVec);
696  break;
697  default:
698  v = u;
699  }
700 }
701 
702 int
703 update_basis(vector<socketstream*> & sock, const VisWinLayout & vwl,
704  Element::Type e, char bType, int bOrder, int mType,
705  Deformation::DefType dType, const DeformationData & defData,
706  bool visualization, int &onlySome)
707 {
708  bool vec = false;
709 
710  Mesh *mesh;
711  ElementMeshStream imesh(e);
712  if (!imesh)
713  {
714  {
715  cerr << "\nProblem with meshstream object\n" << endl;
716  }
717  return 2;
718  }
719  mesh = new Mesh(imesh, 1, 1);
720  int dim = mesh->Dimension();
721 
722  if ( dType != Deformation::INVALID )
723  {
724  Deformation defCoef(dim, dType, defData);
725  mesh->Transform(defCoef);
726  }
727 
728  FiniteElementCollection * FEC = NULL;
729  switch (bType)
730  {
731  case 'h':
732  FEC = new H1_FECollection(bOrder, dim);
733  vec = false;
734  break;
735  case 'p':
736  FEC = new H1Pos_FECollection(bOrder, dim);
737  vec = false;
738  break;
739  case 's':
740  if (bOrder == 1)
741  {
742  FEC = new H1_FECollection(bOrder, dim);
743  }
744  else
745  {
746  FEC = new H1Ser_FECollection(bOrder, dim);
747  }
748  vec = false;
749  break;
750  case 'n':
751  FEC = new ND_FECollection(bOrder, dim);
752  vec = true;
753  break;
754  case 'r':
755  FEC = new RT_FECollection(bOrder-1, dim);
756  vec = true;
757  break;
758  case 'l':
759  FEC = new L2_FECollection(bOrder, dim, BasisType::GaussLegendre,
760  mType);
761  vec = false;
762  break;
763  case 'c':
764  FEC = new CrouzeixRaviartFECollection();
765  break;
766  case 'f':
767  if ( bOrder == 1 )
768  {
769  FEC = new LinearFECollection();
770  }
771  else if ( bOrder == 2 )
772  {
773  FEC = new QuadraticFECollection();
774  }
775  else if ( bOrder == 3 )
776  {
777  FEC = new CubicFECollection();
778  }
779  break;
780  case 'g':
781  if ( bOrder == 1 )
782  {
784  }
785  else if ( bOrder == 2 )
786  {
788  }
789  break;
790  }
791  if ( FEC == NULL)
792  {
793  delete mesh;
794  return 1;
795  }
796 
797  FiniteElementSpace FESpace(mesh, FEC);
798 
799  int ndof = FESpace.GetVSize();
800 
801  Array<int> vdofs;
802  FESpace.GetElementVDofs(0,vdofs);
803 
804  char vishost[] = "localhost";
805  int visport = 19916;
806 
807  int offx = vwl.w+10, offy = vwl.h+45; // window offsets
808 
809  for (unsigned int i=0; i<sock.size(); i++)
810  {
811  *sock[i] << "keys q";
812  delete sock[i];
813  }
814 
815  sock.resize(ndof);
816  for (int i=0; i<ndof; i++)
817  {
818  sock[i] = new socketstream; sock[i]->precision(8);
819  }
820 
821  GridFunction ** x = new GridFunction*[ndof];
822  for (int i=0; i<ndof; i++)
823  {
824  x[i] = new GridFunction(&FESpace);
825  *x[i] = 0.0;
826  if ( vdofs[i] < 0 )
827  {
828  (*x[i])(-1-vdofs[i]) = -1.0;
829  }
830  else
831  {
832  (*x[i])(vdofs[i]) = 1.0;
833  }
834  }
835 
836  int ref = 0;
837  int exOrder = 0;
838  if ( bType == 'n' ) { exOrder++; }
839  if ( bType == 'r' ) { exOrder += 2; }
840  while ( 1<<ref < bOrder + exOrder || ref == 0 )
841  {
842  mesh->UniformRefinement();
843  FESpace.Update();
844 
845  for (int i=0; i<ndof; i++)
846  {
847  x[i]->Update();
848  }
849  ref++;
850  }
851 
852  int stopAt = ndof;
853  if (ndof > 25 && onlySome == -1)
854  {
855  cout << endl;
856  cout << "There are more than 25 windows to open.\n"
857  << "Only showing Dofs 1-10 to avoid crashing.\n"
858  << "Use the option -only N to show Dofs N to N+9 instead.\n";
859  onlySome = 1;
860  }
861  for (int i = 0; i < stopAt; i++)
862  {
863  if (i ==0 && onlySome > 0 && onlySome <ndof)
864  {
865  i = onlySome-1;
866  stopAt = min(ndof,onlySome+9);
867  }
868 
869  ostringstream oss;
870  oss << "DoF " << i + 1;
871  if (visualization)
872  {
873  VisualizeField(*sock[i], vishost, visport, *x[i], oss.str().c_str(),
874  (i % vwl.nx) * offx, ((i / vwl.nx) % vwl.ny) * offy,
875  vwl.w, vwl.h,
876  "aaAc", vec);
877  }
878  }
879 
880  for (int i=0; i<ndof; i++)
881  {
882  delete x[i];
883  }
884  delete [] x;
885 
886  delete FEC;
887  delete mesh;
888 
889  return 0;
890 }
bool elemIs3D(const Element::Type &eType)
Version of QuadraticDiscont2DFECollection with dofs in the Gaussian points.
Definition: fe_coll.hpp:1058
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
virtual void Update(bool want_transform=true)
Reflect changes in the mesh: update number of DOFs, etc. Also, calculate GridFunction transformation ...
Definition: fespace.cpp:3336
Base class for vector Coefficients that optionally depend on time and space.
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
int Dimension() const
Definition: mesh.hpp:1047
bool elemIs1D(const Element::Type &eType)
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
int offx
Piecewise-(bi/tri)linear continuous finite elements.
Definition: fe_coll.hpp:695
void Transform(void(*f)(const Vector &, Vector &))
Definition: mesh.cpp:11691
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150
STL namespace.
Piecewise-(bi)cubic continuous finite elements.
Definition: fe_coll.hpp:775
string mapTypeStr(int mType)
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
constexpr char vishost[]
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9878
constexpr int visport
bool basisIs1D(char bType)
virtual void Eval(DenseMatrix &M, ElementTransformation &T, const IntegrationRule &ir)
Evaluate the vector coefficient in the element described by T at all points of ir, storing the result in M.
Type
Constants for the classes derived from Element.
Definition: element.hpp:41
string elemTypeStr(const Element::Type &eType)
Crouzeix-Raviart nonconforming elements in 2D.
Definition: fe_coll.hpp:806
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition: fe_coll.hpp:373
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: fespace.cpp:281
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
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
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition: gridfunc.cpp:165
bool basisIs3D(char bType)
Version of LinearDiscont2DFECollection with dofs in the Gaussian points.
Definition: fe_coll.hpp:976
Arbitrary order H1-conforming (continuous) finite elements with positive basis functions.
Definition: fe_coll.hpp:292
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:253
Class for integration point with weight.
Definition: intrules.hpp:25
int dim
Definition: ex24.cpp:53
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:587
int main(int argc, char *argv[])
bool elemIs2D(const Element::Type &eType)
Piecewise-(bi)quadratic continuous finite elements.
Definition: fe_coll.hpp:723
bool basisIs2D(char bType)
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition: fe_coll.hpp:447
Vector data type.
Definition: vector.hpp:60
int offy
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:252
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
int update_basis(vector< socketstream *> &sock, const VisWinLayout &vwl, Element::Type e, char bType, int bOrder, int mType, Deformation::DefType dType, const DeformationData &defData, bool visualization, int &onlySome)
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:320
string basisTypeStr(char bType)