41 static int order_ = 3;
44 static double a_ = 1.0;
45 static double b_ = 1.0;
46 static double c_ = 3.0;
48 void pts(
int iphi,
int t,
double x[]);
51 int main(
int argc,
char *argv[])
53 int ser_ref_levels = 0;
57 bool visualization =
true;
60 args.
AddOption(&nz_,
"-nz",
"--num-elements-z",
61 "Number of elements in z-direction.");
62 args.
AddOption(&nt_,
"-nt",
"--num-twists",
63 "Number of node positions to twist the top of the mesh.");
64 args.
AddOption(&order_,
"-o",
"--mesh-order",
65 "Order (polynomial degree) of the mesh elements.");
66 args.
AddOption(&ser_ref_levels,
"-rs",
"--refine-serial",
67 "Number of times to refine the mesh uniformly in serial.");
69 "Width of the base in x-direction.");
71 "Width of the base in y-direction.");
73 "Height in z-direction.");
74 args.
AddOption(&el_type,
"-e",
"--element-type",
75 "Element type: 4 - Tetrahedron, 6 - Wedge, 8 - Hexahedron.");
76 args.
AddOption(&per_mesh,
"-pm",
"--periodic-mesh",
77 "-no-pm",
"--non-periodic-mesh",
78 "Enforce periodicity in z-direction "
79 "(requires discontinuous mesh).");
80 args.
AddOption(&dg_mesh,
"-dm",
"--discont-mesh",
"-cm",
"--cont-mesh",
81 "Use discontinuous or continuous space for the mesh nodes.");
82 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
84 "Enable or disable GLVis visualization.");
97 el_type_ = Element::TETRAHEDRON;
100 el_type_ = Element::WEDGE;
103 el_type_ = Element::HEXAHEDRON;
106 cout <<
"Unsupported element type" << endl;
112 Mesh *mesh =
new Mesh(1, 1, nz_, el_type_,
false, a_, b_, c_,
false);
115 if (order_ > 1 || dg_mesh || per_mesh)
117 mesh->
SetCurvature(order_, dg_mesh || per_mesh, 3, Ordering::byVDIM);
127 if (nt_ % 2 == 1 && fabs(a_ - b_) > 1e-6 * a_)
129 cout <<
"Base is rectangular so number of shifts must be even "
130 <<
"for a periodic mesh!" << endl;
135 if (nt_ % 2 == 1 && (el_type_ == Element::TETRAHEDRON ||
136 el_type_ == Element::WEDGE))
138 cout <<
"Diagonal cuts on the base and top must line up "
139 <<
"for a periodic mesh!" << endl;
144 int noff = (nt_ >= 0) ? 0 : (nnode * (1 - nt_ / nnode));
147 for (
int i = 0; i < v2v.Size() - nnode; i++)
152 switch ((noff + nt_) % nnode)
155 v2v[v2v.Size() - nnode + 0] = 0;
156 v2v[v2v.Size() - nnode + 1] = 1;
157 v2v[v2v.Size() - nnode + 2] = 2;
158 v2v[v2v.Size() - nnode + 3] = 3;
161 v2v[v2v.Size() - nnode + 0] = 2;
162 v2v[v2v.Size() - nnode + 1] = 0;
163 v2v[v2v.Size() - nnode + 2] = 3;
164 v2v[v2v.Size() - nnode + 3] = 1;
167 v2v[v2v.Size() - nnode + 0] = 3;
168 v2v[v2v.Size() - nnode + 1] = 2;
169 v2v[v2v.Size() - nnode + 2] = 1;
170 v2v[v2v.Size() - nnode + 3] = 0;
173 v2v[v2v.Size() - nnode + 0] = 1;
174 v2v[v2v.Size() - nnode + 1] = 3;
175 v2v[v2v.Size() - nnode + 2] = 0;
176 v2v[v2v.Size() - nnode + 3] = 2;
180 for (
int i = 0; i < mesh->
GetNE(); i++)
185 for (
int j = 0; j < nv; j++)
191 for (
int i = 0; i < mesh->
GetNBE(); i++)
196 for (
int j = 0; j < nv; j++)
208 for (
int lev = 0; lev < ser_ref_levels; lev++)
216 if (el_type_ == Element::TETRAHEDRON)
220 else if (el_type_ == Element::WEDGE)
222 oss <<
"twist-wedge";
228 oss <<
"-o" << order_ <<
"-s" << nt_;
229 if (ser_ref_levels > 0)
231 oss <<
"-r" << ser_ref_levels;
247 ofstream ofs(oss.str().c_str());
259 sol_sock.precision(8);
260 sol_sock <<
"mesh\n" << *mesh << flush;
271 double phi = 0.5 * M_PI * nt_ * z / c_;
272 double cp = cos(phi);
273 double sp = sin(phi);
275 p[0] = 0.5 * a_ + (x[0] - 0.5 * a_) * cp - (x[1] - 0.5 * b_) * sp;
276 p[1] = 0.5 * b_ + (x[0] - 0.5 * a_) * sp + (x[1] - 0.5 * b_) * cp;
virtual void Print(std::ostream &out=mfem::out) const
void trans(const Vector &u, Vector &x)
virtual void GetVertices(Array< int > &v) const =0
Returns element's vertices.
int GetNBE() const
Returns number of boundary elements.
void Transform(void(*f)(const Vector &, Vector &))
int GetNE() const
Returns number of elements.
int main(int argc, char *argv[])
void RemoveInternalBoundaries()
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Type
Constants for the classes derived from Element.
const Element * GetElement(int i) const
void PrintUsage(std::ostream &out) const
Print the usage message.
double p(const Vector &x, double t)
void RemoveUnusedVertices()
Remove unused vertices and rebuild mesh connectivity.
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...
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
void PrintOptions(std::ostream &out) const
Print the options.
virtual int GetNVertices() const =0
void pts(int iphi, int t, double x[])
Abstract data type element.
const Element * GetBdrElement(int i) const
bool Good() const
Return true if the command line options were parsed successfully.