MFEM  v4.1.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
display-basis.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, 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)");
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 < 6 )
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  }
276  cout << "enter new element type --> " << flush;
277  cin >> eInt;
278  if ( eInt <= 0 || eInt > 5 )
279  {
280  cout << "invalid element type \"" << eInt << "\"" << endl << flush;
281  }
282  else if ( (elemIs1D((Element::Type)eInt) && basisIs1D(bType)) ||
283  (elemIs2D((Element::Type)eInt) && basisIs2D(bType)) ||
284  (elemIs3D((Element::Type)eInt) && basisIs3D(bType)) )
285  {
286  if ( (elemIs1D((Element::Type)eInt) && !elemIs1D(eType)) ||
287  (elemIs2D((Element::Type)eInt) && !elemIs2D(eType)) ||
288  (elemIs3D((Element::Type)eInt) && !elemIs3D(eType)) )
289  {
290  dType = Deformation::INVALID;
291  }
292  eType = (Element::Type)eInt;
293 
294  print_char = true;
295  }
296  else
297  {
298  cout << "invalid element type \"" << eInt <<
299  "\" for basis type \"" << basisTypeStr(bType) << "\"." << endl;
300  }
301  }
302  if (mk == 'b')
303  {
304  char bChar = 0;
305  cout << "valid basis types:\n";
306  cout << "h) H1 Finite Element\n";
307  cout << "p) H1 Positive Finite Element\n";
308  if ( elemIs2D(eType) || elemIs3D(eType) )
309  {
310  cout << "s) H1 Serendipity Finite Element\n";
311  cout << "n) Nedelec Finite Element\n";
312  cout << "r) Raviart-Thomas Finite Element\n";
313  }
314  cout << "l) L2 Finite Element\n";
315  if ( elemIs1D(eType) || elemIs2D(eType) )
316  {
317  cout << "c) Crouzeix-Raviart Finite Element\n";
318  }
319  cout << "f) Fixed Order Continuous Finite Element\n";
320  if ( elemIs2D(eType) )
321  {
322  cout << "g) Gauss Discontinuous Finite Element\n";
323  }
324  cout << "enter new basis type --> " << flush;
325  cin >> bChar;
326  if (bChar == 'h' || bChar == 'p' || bChar == 'l' || bChar == 'f' ||
327  bChar == 's' ||
328  ((bChar == 'n' || bChar == 'r') && (elemIs2D(eType) || elemIs3D(eType))) ||
329  (bChar == 'c' && (elemIs1D(eType) || elemIs2D(eType))) ||
330  (bChar == 'g' && elemIs2D(eType)))
331  {
332  bType = bChar;
333  if ( bType == 'h' )
334  {
335  mType = FiniteElement::VALUE;
336  }
337  else if ( bType == 'p' )
338  {
339  mType = FiniteElement::VALUE;
340  }
341  else if (bType == 's')
342  {
343  mType = FiniteElement::VALUE;
344  }
345  else if ( bType == 'n' )
346  {
347  mType = FiniteElement::H_CURL;
348  }
349  else if ( bType == 'r' )
350  {
351  mType = FiniteElement::H_DIV;
352  }
353  else if ( bType == 'l' )
354  {
355  if ( mType != FiniteElement::VALUE &&
356  mType != FiniteElement::INTEGRAL )
357  {
358  mType = FiniteElement::VALUE;
359  }
360  }
361  else if ( bType == 'c' )
362  {
363  bOrder = 1;
364  mType = FiniteElement::VALUE;
365  }
366  else if ( bType == 'f' )
367  {
368  if ( bOrder < 1 || bOrder > 3)
369  {
370  bOrder = 1;
371  }
372  mType = FiniteElement::VALUE;
373  }
374  else if ( bType == 'g' )
375  {
376  if ( bOrder < 1 || bOrder > 2)
377  {
378  bOrder = 1;
379  }
380  mType = FiniteElement::VALUE;
381  }
382  print_char = true;
383  }
384  else
385  {
386  cout << "invalid basis type \"" << bChar << "\"." << endl;
387  }
388  }
389  if (mk == 'm' && bType == 'l')
390  {
391  int mInt = 0;
392  cout << "valid map types:\n"
393  "0) VALUE\n"
394  "1) INTEGRAL\n";
395  cout << "enter new map type --> " << flush;
396  cin >> mInt;
397  if (mInt >=0 && mInt <= 1)
398  {
399  mType = mInt;
400  print_char = true;
401  }
402  else
403  {
404  cout << "invalid map type \"" << mInt << "\"." << endl;
405  }
406  }
407  if (mk == 'o')
408  {
409  int oInt = 1;
410  int oMin = ( bType == 'h' || bType == 'p' || bType == 'n' ||
411  bType == 'f' || bType == 'g' || bType == 's')?1:0;
412  int oMax = -1;
413  switch (bType)
414  {
415  case 'g':
416  oMax = 2;
417  break;
418  case 'f':
419  oMax = 3;
420  break;
421  default:
422  oMax = -1;
423  }
424  cout << "basis function order must be >= " << oMin;
425  if ( oMax >= 0 )
426  {
427  cout << " and <= " << oMax;
428  }
429  cout << endl;
430  cout << "enter new basis function order --> " << flush;
431  cin >> oInt;
432  if ( oInt >= oMin && oInt <= (oMax>=0)?oMax:oInt )
433  {
434  bOrder = oInt;
435  print_char = true;
436  }
437  else
438  {
439  cout << "invalid basis order \"" << oInt << "\"." << endl;
440  }
441  }
442  if (mk == 't')
443  {
444  cout << "transformation options:\n";
445  cout << "r) reset to reference element\n";
446  cout << "u) uniform scaling\n";
447  if ( elemIs2D(eType) || elemIs3D(eType) )
448  {
449  cout << "c) compression\n";
450  cout << "s) shear\n";
451  }
452  cout << "enter transformation type --> " << flush;
453  char tk;
454  cin >> tk;
455  if (tk == 'r')
456  {
457  dType = Deformation::INVALID;
458  }
459  else if (tk == 'u')
460  {
461  cout << "enter scaling constant --> " << flush;
462  cin >> defData.uniformScale;
463  if ( defData.uniformScale > 0.0 )
464  {
465  dType = Deformation::UNIFORM;
466  }
467  }
468  else if (tk == 'c' && !elemIs1D(eType))
469  {
470  int dim = elemIs2D(eType)?2:3;
471  cout << "enter compression factor --> " << flush;
472  cin >> defData.squeezeFactor;
473  cout << "enter compression axis (0-" << dim-1 << ") --> " << flush;
474  cin >> defData.squeezeAxis;
475 
476  if ( defData.squeezeFactor > 0.0 &&
477  (defData.squeezeAxis >= 0 && defData.squeezeAxis < dim))
478  {
479  dType = Deformation::SQUEEZE;
480  }
481  }
482  else if (tk == 's' && !elemIs1D(eType))
483  {
484  int dim = elemIs2D(eType)?2:3;
485  cout << "enter shear vector (components separated by spaces) --> "
486  << flush;
487  defData.shearVec.SetSize(dim);
488  for (int i=0; i<dim; i++)
489  {
490  cin >> defData.shearVec[i];
491  }
492  cout << "enter shear axis (0-" << dim-1 << ") --> " << flush;
493  cin >> defData.shearAxis;
494 
495  if ( defData.shearAxis >= 0 && defData.shearAxis < dim )
496  {
497  dType = Deformation::SHEAR;
498  }
499  }
500  }
501  }
502 
503  // Cleanup
504  for (unsigned int i=0; i<sock.size(); i++)
505  {
506  delete sock[i];
507  }
508 
509  // Exit
510  return 0;
511 }
512 
513 string elemTypeStr(const Element::Type & eType)
514 {
515  switch (eType)
516  {
517  case Element::POINT:
518  return "POINT";
519  case Element::SEGMENT:
520  return "SEGMENT";
521  case Element::TRIANGLE:
522  return "TRIANGLE";
523  case Element::QUADRILATERAL:
524  return "QUADRILATERAL";
525  case Element::TETRAHEDRON:
526  return "TETRAHEDRON";
527  case Element::HEXAHEDRON:
528  return "HEXAHEDRON";
529  default:
530  return "INVALID";
531  };
532 }
533 
534 bool
535 elemIs1D(const Element::Type & eType)
536 {
537  return eType == Element::SEGMENT;
538 }
539 
540 bool
541 elemIs2D(const Element::Type & eType)
542 {
543  return eType == Element::TRIANGLE || eType == Element::QUADRILATERAL;
544 }
545 
546 bool
547 elemIs3D(const Element::Type & eType)
548 {
549  return eType == Element::TETRAHEDRON || eType == Element::HEXAHEDRON;
550 }
551 
552 string
553 basisTypeStr(char bType)
554 {
555  switch (bType)
556  {
557  case 'h':
558  return "Continuous (H1)";
559  case 'p':
560  return "Continuous Positive (H1)";
561  case 's':
562  return "Continuous Serendipity (H1)";
563  case 'n':
564  return "Nedelec";
565  case 'r':
566  return "Raviart-Thomas";
567  case 'l':
568  return "Discontinuous (L2)";
569  case 'f':
570  return "Fixed Order Continuous";
571  case 'g':
572  return "Gaussian Discontinuous";
573  case 'c':
574  return "Crouzeix-Raviart";
575  default:
576  return "INVALID";
577  };
578 }
579 
580 bool
581 basisIs1D(char bType)
582 {
583  return bType == 'h' || bType == 'p' || bType == 'l' || bType == 'c' ||
584  bType == 'f';
585 }
586 
587 bool
588 basisIs2D(char bType)
589 {
590  return bType == 'h' || bType == 'p' || bType == 'n' || bType == 'r' ||
591  bType == 'l' || bType == 'c' || bType == 'f' || bType == 'g' ||
592  bType == 's';
593 }
594 
595 bool
596 basisIs3D(char bType)
597 {
598  return bType == 'h' || bType == 'p' || bType == 'n' || bType == 'r' ||
599  bType == 'f' || bType == 'l';
600 }
601 
602 string
603 mapTypeStr(int mType)
604 {
605  switch (mType)
606  {
607  case FiniteElement::VALUE:
608  return "VALUE";
609  case FiniteElement::H_CURL:
610  return "H_CURL";
611  case FiniteElement::H_DIV:
612  return "H_DIV";
613  case FiniteElement::INTEGRAL:
614  return "INTEGRAL";
615  default:
616  return "INVALID";
617  }
618 }
619 
620 void
622  const IntegrationPoint &ip)
623 {
624  Vector u(dim_);
625  T.Transform(ip, u);
626 
627  switch (dim_)
628  {
629  case 1:
630  Def1D(u, v);
631  break;
632  case 2:
633  Def2D(u, v);
634  break;
635  case 3:
636  Def3D(u, v);
637  break;
638  }
639 }
640 
641 void
642 Deformation::Def1D(const Vector & u, Vector & v)
643 {
644  v = u;
645  if ( dType_ == UNIFORM )
646  {
647  v *= data_.uniformScale;
648  }
649 }
650 
651 void
652 Deformation::Def2D(const Vector & u, Vector & v)
653 {
654  switch (dType_)
655  {
656  case UNIFORM:
657  v = u;
658  v *= data_.uniformScale;
659  break;
660  case SQUEEZE:
661  v = u;
662  v[ data_.squeezeAxis ] /= data_.squeezeFactor;
663  v[(data_.squeezeAxis+1)%2] *= data_.squeezeFactor;
664  break;
665  case SHEAR:
666  v = u;
667  v.Add(v[data_.shearAxis], data_.shearVec);
668  break;
669  default:
670  v = u;
671  }
672 }
673 
674 void
675 Deformation::Def3D(const Vector & u, Vector & v)
676 {
677  switch (dType_)
678  {
679  case UNIFORM:
680  v = u;
681  v *= data_.uniformScale;
682  break;
683  case SQUEEZE:
684  v = u;
685  v[ data_.squeezeAxis ] /= data_.squeezeFactor;
686  v[(data_.squeezeAxis+1)%2] *= sqrt(data_.squeezeFactor);
687  v[(data_.squeezeAxis+2)%2] *= sqrt(data_.squeezeFactor);
688  break;
689  case SHEAR:
690  v = u;
691  v.Add(v[data_.shearAxis], data_.shearVec);
692  break;
693  default:
694  v = u;
695  }
696 }
697 
698 int
699 update_basis(vector<socketstream*> & sock, const VisWinLayout & vwl,
700  Element::Type e, char bType, int bOrder, int mType,
701  Deformation::DefType dType, const DeformationData & defData,
702  bool visualization, int &onlySome)
703 {
704  bool vec = false;
705 
706  Mesh *mesh;
707  ElementMeshStream imesh(e);
708  if (!imesh)
709  {
710  {
711  cerr << "\nProblem with meshstream object\n" << endl;
712  }
713  return 2;
714  }
715  mesh = new Mesh(imesh, 1, 1);
716  int dim = mesh->Dimension();
717 
718  if ( dType != Deformation::INVALID )
719  {
720  Deformation defCoef(dim, dType, defData);
721  mesh->Transform(defCoef);
722  }
723 
724  FiniteElementCollection * FEC = NULL;
725  switch (bType)
726  {
727  case 'h':
728  FEC = new H1_FECollection(bOrder, dim);
729  vec = false;
730  break;
731  case 'p':
732  FEC = new H1Pos_FECollection(bOrder, dim);
733  vec = false;
734  break;
735  case 's':
736  if (bOrder == 1)
737  {
738  FEC = new H1_FECollection(bOrder, dim);
739  }
740  else
741  {
742  FEC = new H1Ser_FECollection(bOrder, dim);
743  }
744  vec = false;
745  break;
746  case 'n':
747  FEC = new ND_FECollection(bOrder, dim);
748  vec = true;
749  break;
750  case 'r':
751  FEC = new RT_FECollection(bOrder-1, dim);
752  vec = true;
753  break;
754  case 'l':
755  FEC = new L2_FECollection(bOrder, dim, BasisType::GaussLegendre,
756  mType);
757  vec = false;
758  break;
759  case 'c':
760  FEC = new CrouzeixRaviartFECollection();
761  break;
762  case 'f':
763  if ( bOrder == 1 )
764  {
765  FEC = new LinearFECollection();
766  }
767  else if ( bOrder == 2 )
768  {
769  FEC = new QuadraticFECollection();
770  }
771  else if ( bOrder == 3 )
772  {
773  FEC = new CubicFECollection();
774  }
775  break;
776  case 'g':
777  if ( bOrder == 1 )
778  {
780  }
781  else if ( bOrder == 2 )
782  {
784  }
785  break;
786  }
787  if ( FEC == NULL)
788  {
789  delete mesh;
790  return 1;
791  }
792 
793  FiniteElementSpace FESpace(mesh, FEC);
794 
795  int ndof = FESpace.GetVSize();
796 
797  Array<int> vdofs;
798  FESpace.GetElementVDofs(0,vdofs);
799 
800  char vishost[] = "localhost";
801  int visport = 19916;
802 
803  int offx = vwl.w+10, offy = vwl.h+45; // window offsets
804 
805  for (unsigned int i=0; i<sock.size(); i++)
806  {
807  *sock[i] << "keys q";
808  delete sock[i];
809  }
810 
811  sock.resize(ndof);
812  for (int i=0; i<ndof; i++)
813  {
814  sock[i] = new socketstream; sock[i]->precision(8);
815  }
816 
817  GridFunction ** x = new GridFunction*[ndof];
818  for (int i=0; i<ndof; i++)
819  {
820  x[i] = new GridFunction(&FESpace);
821  *x[i] = 0.0;
822  if ( vdofs[i] < 0 )
823  {
824  (*x[i])(-1-vdofs[i]) = -1.0;
825  }
826  else
827  {
828  (*x[i])(vdofs[i]) = 1.0;
829  }
830  }
831 
832  int ref = 0;
833  int exOrder = 0;
834  if ( bType == 'n' ) { exOrder++; }
835  if ( bType == 'r' ) { exOrder += 2; }
836  while ( 1<<ref < bOrder + exOrder || ref == 0 )
837  {
838  mesh->UniformRefinement();
839  FESpace.Update();
840 
841  for (int i=0; i<ndof; i++)
842  {
843  x[i]->Update();
844  }
845  ref++;
846  }
847 
848  int stopAt = ndof;
849  if (ndof > 25 && onlySome == -1)
850  {
851  cout << endl;
852  cout << "There are more than 25 windows to open.\n"
853  << "Only showing Dofs 1-10 to avoid crashing.\n"
854  << "Use the option -only N to show Dofs N to N+9 instead.\n";
855  onlySome = 1;
856  }
857  for (int i = 0; i < stopAt; i++)
858  {
859  if (i ==0 && onlySome > 0 && onlySome <ndof)
860  {
861  i = onlySome-1;
862  stopAt = min(ndof,onlySome+9);
863  }
864 
865  ostringstream oss;
866  oss << "DoF " << i + 1;
867  if (visualization)
868  {
869  VisualizeField(*sock[i], vishost, visport, *x[i], oss.str().c_str(),
870  (i % vwl.nx) * offx, ((i / vwl.nx) % vwl.ny) * offy,
871  vwl.w, vwl.h,
872  "aaAc", vec);
873  }
874  }
875 
876  for (int i=0; i<ndof; i++)
877  {
878  delete x[i];
879  }
880  delete [] x;
881 
882  delete FEC;
883  delete mesh;
884 
885  return 0;
886 }
bool elemIs3D(const Element::Type &eType)
Version of QuadraticDiscont2DFECollection with dofs in the Gaussian points.
Definition: fe_coll.hpp:672
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:381
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:27
virtual void Update(bool want_transform=true)
Definition: fespace.cpp:2057
bool elemIs1D(const Element::Type &eType)
void GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: fespace.cpp:172
int offx
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)
Piecewise-(bi)linear continuous finite elements.
Definition: fe_coll.hpp:343
void Transform(void(*f)(const Vector &, Vector &))
Definition: mesh.cpp:9864
int main(int argc, char *argv[])
Definition: ex1.cpp:62
Piecewise-(bi)cubic continuous finite elements.
Definition: fe_coll.hpp:415
string mapTypeStr(int mType)
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:7982
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.
Definition: coefficient.cpp:90
Type
Constants for the classes derived from Element.
Definition: element.hpp:41
string elemTypeStr(const Element::Type &eType)
int Dimension() const
Definition: mesh.hpp:744
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
Crouzeix-Raviart nonconforming elements in 2D.
Definition: fe_coll.hpp:442
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition: fe_coll.hpp:191
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
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:76
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition: gridfunc.cpp:152
bool basisIs3D(char bType)
Version of LinearDiscont2DFECollection with dofs in the Gaussian points.
Definition: fe_coll.hpp:596
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:190
Class for integration point with weight.
Definition: intrules.hpp:25
int dim
Definition: ex24.cpp:43
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
bool elemIs2D(const Element::Type &eType)
Piecewise-(bi)quadratic continuous finite elements.
Definition: fe_coll.hpp:368
bool basisIs2D(char bType)
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition: fe_coll.hpp:251
Vector data type.
Definition: vector.hpp:48
void VisualizeField(socketstream &sock, const char *vishost, int visport, GridFunction &gf, const char *title, int x, int y, int w, int h, const char *keys, bool vec)
Definition: fem_extras.cpp:94
virtual void Transform(const IntegrationPoint &, Vector &)=0
int offy
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:83
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:143
string basisTypeStr(char bType)
bool Good() const
Definition: optparser.hpp:125