MFEM  v4.3.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
twist.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2021, 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 // Twist Miniapp: Generate simple twisted periodic meshes
14 // -------------------------------------------------------
15 //
16 // This miniapp generates simple periodic meshes to demonstrate MFEM's handling
17 // of periodic domains. MFEM's strategy is to use a discontinuous vector field
18 // to define the mesh coordinates on a topologically periodic mesh. It works by
19 // defining a stack of individual elements and stitching together the top and
20 // bottom of the mesh. The stack can also be twisted so that the vertices of the
21 // bottom and top can be joined with any integer offset (for tetrahedral and
22 // wedge meshes only even offsets are supported).
23 //
24 // Compile with: make twist
25 //
26 // Sample runs: twist
27 // twist -no-pm
28 // twist -nt -2 -no-pm
29 // twist -nt 2 -e 4
30 // twist -nt 2 -e 6
31 // twist -nt 3 -e 8
32 //
33 #include "mfem.hpp"
34 #include <fstream>
35 #include <iostream>
36 
37 using namespace std;
38 using namespace mfem;
39 
40 static Element::Type el_type_ = Element::WEDGE;
41 static int order_ = 3;
42 static int nz_ = 3;
43 static int nt_ = 2;
44 static double a_ = 1.0;
45 static double b_ = 1.0;
46 static double c_ = 3.0;
47 
48 void pts(int iphi, int t, double x[]);
49 void trans(const Vector &x, Vector &p);
50 
51 int main(int argc, char *argv[])
52 {
53  int ser_ref_levels = 0;
54  int el_type = 8;
55  bool per_mesh = true;
56  bool dg_mesh = false;
57  bool visualization = true;
58 
59  OptionsParser args(argc, argv);
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.");
68  args.AddOption(&a_, "-a", "--base-x",
69  "Width of the base in x-direction.");
70  args.AddOption(&b_, "-b", "--base-y",
71  "Width of the base in y-direction.");
72  args.AddOption(&c_, "-c", "--height",
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",
83  "--no-visualization",
84  "Enable or disable GLVis visualization.");
85  args.Parse();
86  if (!args.Good())
87  {
88  args.PrintUsage(cout);
89  return 1;
90  }
91  args.PrintOptions(cout);
92 
93  // The output mesh could be tetrahedra, hexahedra, or prisms
94  switch (el_type)
95  {
96  case 4:
97  el_type_ = Element::TETRAHEDRON;
98  break;
99  case 6:
100  el_type_ = Element::WEDGE;
101  break;
102  case 8:
103  el_type_ = Element::HEXAHEDRON;
104  break;
105  default:
106  cout << "Unsupported element type" << endl;
107  exit(1);
108  break;
109  }
110 
111  // Define the mesh
112  Mesh mesh = Mesh::MakeCartesian3D(1, 1, nz_, el_type_, a_, b_, c_, false);
113 
114  // Promote to high order mesh and transform into a twisted shape
115  if (order_ > 1 || dg_mesh || per_mesh)
116  {
117  mesh.SetCurvature(order_, dg_mesh || per_mesh, 3, Ordering::byVDIM);
118  }
119  if (nt_ != 0 )
120  {
121  mesh.Transform(trans);
122  }
123 
124  while (per_mesh)
125  {
126  // Verify geometric compatibility
127  if (nt_ % 2 == 1 && fabs(a_ - b_) > 1e-6 * a_)
128  {
129  cout << "Base is rectangular so number of shifts must be even "
130  << "for a periodic mesh!" << endl;
131  exit(1);
132  }
133 
134  // Verify topological compatibility
135  if (nt_ % 2 == 1 && (el_type_ == Element::TETRAHEDRON ||
136  el_type_ == Element::WEDGE))
137  {
138  cout << "Diagonal cuts on the base and top must line up "
139  << "for a periodic mesh!" << endl;
140  exit(1);
141  }
142 
143  int nnode = 4;
144  int noff = (nt_ >= 0) ? 0 : (nnode * (1 - nt_ / nnode));
145 
146  Array<int> v2v(mesh.GetNV());
147  for (int i = 0; i < v2v.Size() - nnode; i++)
148  {
149  v2v[i] = i;
150  }
151  // identify vertices at the extremes of the stack
152  switch ((noff + nt_) % nnode)
153  {
154  case 0:
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;
159  break;
160  case 1:
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;
165  break;
166  case 2:
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;
171  break;
172  case 3:
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;
177  break;
178  }
179  // renumber elements
180  for (int i = 0; i < mesh.GetNE(); i++)
181  {
182  Element *el = mesh.GetElement(i);
183  int *v = el->GetVertices();
184  int nv = el->GetNVertices();
185  for (int j = 0; j < nv; j++)
186  {
187  v[j] = v2v[v[j]];
188  }
189  }
190  // renumber boundary elements
191  for (int i = 0; i < mesh.GetNBE(); i++)
192  {
193  Element *el = mesh.GetBdrElement(i);
194  int *v = el->GetVertices();
195  int nv = el->GetNVertices();
196  for (int j = 0; j < nv; j++)
197  {
198  v[j] = v2v[v[j]];
199  }
200  }
201  mesh.RemoveUnusedVertices();
203 
204  break;
205  }
206 
207  // Refine the mesh if desired
208  for (int lev = 0; lev < ser_ref_levels; lev++)
209  {
210  mesh.UniformRefinement();
211  }
212 
213  // Output the resulting mesh to a file
214  {
215  ostringstream oss;
216  if (el_type_ == Element::TETRAHEDRON)
217  {
218  oss << "twist-tet";
219  }
220  else if (el_type_ == Element::WEDGE)
221  {
222  oss << "twist-wedge";
223  }
224  else
225  {
226  oss << "twist-hex";
227  }
228  oss << "-o" << order_ << "-s" << nt_;
229  if (ser_ref_levels > 0)
230  {
231  oss << "-r" << ser_ref_levels;
232  }
233  if (per_mesh)
234  {
235  oss << "-p";
236  }
237  else if (dg_mesh)
238  {
239  oss << "-d";
240  }
241  else
242  {
243  oss << "-c";
244  }
245 
246  oss << ".mesh";
247  ofstream ofs(oss.str().c_str());
248  ofs.precision(8);
249  mesh.Print(ofs);
250  ofs.close();
251  }
252 
253  // Output the resulting mesh to GLVis
254  if (visualization)
255  {
256  char vishost[] = "localhost";
257  int visport = 19916;
258  socketstream sol_sock(vishost, visport);
259  sol_sock.precision(8);
260  sol_sock << "mesh\n" << mesh << flush;
261  }
262 
263  // Clean up and exit
264  return 0;
265 }
266 
267 void trans(const Vector &x, Vector &p)
268 {
269  double z = x[2];
270  double phi = 0.5 * M_PI * nt_ * z / c_;
271  double cp = cos(phi);
272  double sp = sin(phi);
273 
274  p[0] = 0.5 * a_ + (x[0] - 0.5 * a_) * cp - (x[1] - 0.5 * b_) * sp;
275  p[1] = 0.5 * b_ + (x[0] - 0.5 * a_) * sp + (x[1] - 0.5 * b_) * cp;
276  p[2] = z;
277 }
virtual void Print(std::ostream &out=mfem::out) const
Definition: mesh.hpp:1387
void trans(const Vector &u, Vector &x)
Definition: ex27.cpp:421
virtual void GetVertices(Array< int > &v) const =0
Returns element&#39;s vertices.
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:849
void Transform(void(*f)(const Vector &, Vector &))
Definition: mesh.cpp:11008
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:846
void RemoveInternalBoundaries()
Definition: mesh.cpp:11164
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:150
constexpr char vishost[]
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9143
constexpr int visport
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:4882
Type
Constants for the classes derived from Element.
Definition: element.hpp:41
const Element * GetElement(int i) const
Definition: mesh.hpp:942
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:457
void RemoveUnusedVertices()
Remove unused vertices and rebuild mesh connectivity.
Definition: mesh.cpp:11057
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
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition: mesh.hpp:843
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:327
virtual int GetNVertices() const =0
RefCoord t[3]
Vector data type.
Definition: vector.hpp:60
void pts(int iphi, int t, double x[])
Abstract data type element.
Definition: element.hpp:28
int main()
const Element * GetBdrElement(int i) const
Definition: mesh.hpp:946
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150