29 #include "../common/fem_extras.hpp"
30 #include "../common/mesh_extras.hpp"
36 using namespace mfem::miniapps;
48 struct DeformationData
76 enum DefType {INVALID, UNIFORM, SQUEEZE, SHEAR};
78 Deformation(
int dim, DefType dType,
const DeformationData & data)
90 const DeformationData & data_;
105 int update_basis(vector<socketstream*> & sock,
const VisWinLayout & vwl,
107 Deformation::DefType dType,
const DeformationData & defData,
110 int main(
int argc,
char *argv[])
127 Deformation::DefType dType = Deformation::INVALID;
128 DeformationData defData;
130 bool visualization =
true;
132 vector<socketstream*> sock;
135 args.
AddOption(&eInt,
"-e",
"--elem-type",
136 "Element Type: (1-Segment, 2-Triangle, 3-Quadrilateral, "
137 "4-Tetrahedron, 5-Hexahedron)");
138 args.
AddOption(&bInt,
"-b",
"--basis-type",
139 "Basis Function Type (0-H1, 1-Nedelec, 2-Raviart-Thomas, "
140 "3-L2, 4-Fixed Order Cont.,\n\t5-Gaussian Discontinuous (2D),"
141 " 6-Crouzeix-Raviart)");
142 args.
AddOption(&bOrder,
"-o",
"--order",
"Basis function order");
143 args.
AddOption(&vwl.nx,
"-nx",
"--num-win-x",
144 "Number of Viz windows in X");
145 args.
AddOption(&vwl.ny,
"-ny",
"--num-win-y",
146 "Number of Viz windows in y");
148 "Width of Viz windows");
150 "Height of Viz windows");
151 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
152 "--no-visualization",
153 "Enable or disable GLVis visualization.");
163 if ( eInt > 0 && eInt < 6 )
195 bool print_char =
true;
201 cout <<
"Element Type: " <<
elemTypeStr(eType) << endl;
203 cout <<
"Basis function order: " << bOrder << endl;
204 cout <<
"Map Type: " <<
mapTypeStr(mType) << endl;
206 if (
update_basis(sock, vwl, eType, bType, bOrder, mType,
207 dType, defData, visualization) )
209 cerr <<
"Invalid combination of basis info (try again)" << endl;
212 if (!visualization) {
break; }
216 cout <<
"What would you like to do?\n"
218 "c) Close Windows and Quit\n"
219 "e) Change Element Type\n"
220 "b) Change Basis Type\n";
221 if ( bType ==
'h' || bType ==
'p' || bType ==
'n' || bType ==
'r' ||
222 bType ==
'l' || bType ==
'f' || bType ==
'g' )
224 cout <<
"o) Change Basis Order\n";
227 if ( bType ==
'l' &&
false )
229 cout <<
"m) Change Map Type\n";
231 cout <<
"t) Transform Element\n";
232 cout <<
"--> " << flush;
242 for (
unsigned int i=0; i<sock.size(); i++)
244 *sock[i] <<
"keys q";
251 cout <<
"valid element types:\n";
261 "3) Quadrilateral\n";
269 cout <<
"enter new element type --> " << flush;
271 if ( eInt <= 0 || eInt > 5 )
273 cout <<
"invalid element type \"" << eInt <<
"\"" << endl << flush;
283 dType = Deformation::INVALID;
291 cout <<
"invalid element type \"" << eInt <<
292 "\" for basis type \"" <<
basisTypeStr(bType) <<
"\"." << endl;
298 cout <<
"valid basis types:\n";
299 cout <<
"h) H1 Finite Element\n";
300 cout <<
"p) H1 Positive Finite Element\n";
303 cout <<
"n) Nedelec Finite Element\n";
304 cout <<
"r) Raviart-Thomas Finite Element\n";
306 cout <<
"l) L2 Finite Element\n";
309 cout <<
"c) Crouzeix-Raviart Finite Element\n";
311 cout <<
"f) Fixed Order Continuous Finite Element\n";
314 cout <<
"g) Gauss Discontinuous Finite Element\n";
316 cout <<
"enter new basis type --> " << flush;
318 if ( bChar ==
'h' || bChar ==
'p' || bChar ==
'l' || bChar ==
'f' ||
319 ((bChar ==
'n' || bChar ==
'r') &&
327 mType = FiniteElement::VALUE;
329 else if ( bType ==
'p' )
331 mType = FiniteElement::VALUE;
333 else if ( bType ==
'n' )
335 mType = FiniteElement::H_CURL;
337 else if ( bType ==
'r' )
339 mType = FiniteElement::H_DIV;
341 else if ( bType ==
'l' )
343 if ( mType != FiniteElement::VALUE &&
344 mType != FiniteElement::INTEGRAL )
346 mType = FiniteElement::VALUE;
349 else if ( bType ==
'c' )
352 mType = FiniteElement::VALUE;
354 else if ( bType ==
'f' )
356 if ( bOrder < 1 || bOrder > 3)
360 mType = FiniteElement::VALUE;
362 else if ( bType ==
'g' )
364 if ( bOrder < 1 || bOrder > 2)
368 mType = FiniteElement::VALUE;
374 cout <<
"invalid basis type \"" << bChar <<
"\"." << endl;
377 if (mk ==
'm' && bType ==
'l')
380 cout <<
"valid map types:\n"
383 cout <<
"enter new map type --> " << flush;
385 if (mInt >=0 && mInt <= 1)
392 cout <<
"invalid map type \"" << mInt <<
"\"." << endl;
398 int oMin = ( bType ==
'h' || bType ==
'p' || bType ==
'n' ||
399 bType ==
'f' || bType ==
'g')?1:0;
412 cout <<
"basis function order must be >= " << oMin;
415 cout <<
" and <= " << oMax;
418 cout <<
"enter new basis function order --> " << flush;
420 if ( oInt >= oMin && oInt <= (oMax>=0)?oMax:oInt )
427 cout <<
"invalid basis order \"" << oInt <<
"\"." << endl;
432 cout <<
"transformation options:\n";
433 cout <<
"r) reset to reference element\n";
434 cout <<
"u) uniform scaling\n";
437 cout <<
"c) compression\n";
438 cout <<
"s) shear\n";
440 cout <<
"enter transformation type --> " << flush;
445 dType = Deformation::INVALID;
449 cout <<
"enter scaling constant --> " << flush;
450 cin >> defData.uniformScale;
451 if ( defData.uniformScale > 0.0 )
453 dType = Deformation::UNIFORM;
456 else if (tk ==
'c' && !
elemIs1D(eType))
459 cout <<
"enter compression factor --> " << flush;
460 cin >> defData.squeezeFactor;
461 cout <<
"enter compression axis (0-" << dim-1 <<
") --> " << flush;
462 cin >> defData.squeezeAxis;
464 if ( defData.squeezeFactor > 0.0 &&
465 (defData.squeezeAxis >= 0 && defData.squeezeAxis < dim))
467 dType = Deformation::SQUEEZE;
470 else if (tk ==
's' && !
elemIs1D(eType))
473 cout <<
"enter shear vector (components separated by spaces) --> "
475 defData.shearVec.SetSize(dim);
476 for (
int i=0; i<
dim; i++)
478 cin >> defData.shearVec[i];
480 cout <<
"enter shear axis (0-" << dim-1 <<
") --> " << flush;
481 cin >> defData.shearAxis;
483 if ( defData.shearAxis >= 0 && defData.shearAxis < dim )
485 dType = Deformation::SHEAR;
492 for (
unsigned int i=0; i<sock.size(); i++)
508 case Element::SEGMENT:
511 case Element::TRIANGLE:
514 case Element::QUADRILATERAL:
515 return "QUADRILATERAL";
517 case Element::TETRAHEDRON:
518 return "TETRAHEDRON";
520 case Element::HEXAHEDRON:
532 return eType == Element::SEGMENT;
538 return eType == Element::TRIANGLE || eType == Element::QUADRILATERAL;
544 return eType == Element::TETRAHEDRON || eType == Element::HEXAHEDRON;
553 return "Continuous (H1)";
556 return "Continuous Positive (H1)";
562 return "Raviart-Thomas";
565 return "Discontinuous (L2)";
568 return "Fixed Order Continuous";
571 return "Gaussian Discontinuous";
574 return "Crouzeix-Raviart";
585 return bType ==
'h' || bType ==
'p' || bType ==
'l' || bType ==
'c' ||
592 return bType ==
'h' || bType ==
'p' || bType ==
'n' || bType ==
'r' ||
593 bType ==
'l' || bType ==
'c' || bType ==
'f' || bType ==
'g';
599 return bType ==
'h' || bType ==
'p' || bType ==
'n' || bType ==
'r' ||
600 bType ==
'f' || bType ==
'l';
608 case FiniteElement::VALUE:
611 case FiniteElement::H_CURL:
614 case FiniteElement::H_DIV:
617 case FiniteElement::INTEGRAL:
651 if ( dType_ == UNIFORM )
653 v *= data_.uniformScale;
664 v *= data_.uniformScale;
668 v[ data_.squeezeAxis ] /= data_.squeezeFactor;
669 v[(data_.squeezeAxis+1)%2] *= data_.squeezeFactor;
673 v.
Add(v[data_.shearAxis], data_.shearVec);
687 v *= data_.uniformScale;
691 v[ data_.squeezeAxis ] /= data_.squeezeFactor;
692 v[(data_.squeezeAxis+1)%2] *= sqrt(data_.squeezeFactor);
693 v[(data_.squeezeAxis+2)%2] *= sqrt(data_.squeezeFactor);
697 v.
Add(v[data_.shearAxis], data_.shearVec);
707 Deformation::DefType dType,
const DeformationData & defData,
717 cerr <<
"\nProblem with meshstream object\n" << endl;
721 mesh =
new Mesh(imesh, 1, 1);
724 if ( dType != Deformation::INVALID )
726 Deformation defCoef(dim, dType, defData);
762 else if ( bOrder == 2 )
766 else if ( bOrder == 3 )
776 else if ( bOrder == 2 )
794 char vishost[] =
"localhost";
797 int offx = vwl.w+10, offy = vwl.h+45;
799 for (
unsigned int i=0; i<sock.size(); i++)
801 *sock[i] <<
"keys q";
806 for (
int i=0; i<ndof; i++)
812 for (
int i=0; i<ndof; i++)
818 (*x[i])(-1-vdofs[i]) = -1.0;
822 (*x[i])(vdofs[i]) = 1.0;
828 if ( bType ==
'n' ) { exOrder++; }
829 if ( bType ==
'r' ) { exOrder += 2; }
830 while ( 1<<ref < bOrder + exOrder || ref == 0 )
835 for (
int i=0; i<ndof; i++)
842 for (
int i=0; i<ndof; i++)
845 oss <<
"DoF " << i + 1;
848 VisualizeField(*sock[i], vishost, visport, *x[i], oss.str().c_str(),
849 (i % vwl.nx) * offx, ((i / vwl.nx) % vwl.ny) * offy,
855 for (
int i=0; i<ndof; i++)
bool elemIs3D(const Element::Type &eType)
Version of QuadraticDiscont2DFECollection with dofs in the Gaussian points.
Class for grid function - Vector with associated FE space.
virtual void Update(bool want_transform=true)
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'th element.
Piecewise-(bi)linear continuous finite elements.
void Transform(void(*f)(const Vector &, Vector &))
int main(int argc, char *argv[])
Piecewise-(bi)cubic continuous finite elements.
string mapTypeStr(int mType)
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)
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
bool basisIs1D(char bType)
void VisualizeField(socketstream &sock, const char *vishost, int visport, GridFunction &gf, const char *title, int x, int y, int w, int h, bool vec)
virtual void Eval(DenseMatrix &M, ElementTransformation &T, const IntegrationRule &ir)
Type
Constants for the classes derived from Element.
string elemTypeStr(const Element::Type &eType)
void PrintUsage(std::ostream &out) const
Crouzeix-Raviart nonconforming elements in 2D.
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
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)
void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
bool basisIs3D(char bType)
Version of LinearDiscont2DFECollection with dofs in the Gaussian points.
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Class for integration point with weight.
void PrintOptions(std::ostream &out) const
bool elemIs2D(const Element::Type &eType)
Piecewise-(bi)quadratic continuous finite elements.
bool basisIs2D(char bType)
Arbitrary order H(curl)-conforming Nedelec finite elements.
Arbitrary order H1-conforming (continuous) finite elements.
Arbitrary order "L2-conforming" discontinuous finite elements.
string basisTypeStr(char bType)