MFEM  v4.6.0
Finite element discretization library
toroid.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 // Toroid Miniapp: Generate simple toroidal meshes
14 // ------------------------------------------------
15 //
16 // This miniapp generates two types of toroidal meshes; one with triangular
17 // cross sections and one with square cross sections. It works by defining a
18 // stack of individual elements and bending them so that the bottom and top of
19 // the stack can be joined to form a torus. The stack can also be twisted so
20 // that the vertices of the bottom and top can be joined with any integer
21 // offset.
22 //
23 // Compile with: make toroid
24 //
25 // Sample runs: toroid
26 // toroid -nphi 6
27 // toroid -ns 1
28 // toroid -ns 0 -t0 -30
29 // toroid -R 2 -r 1 -ns 3
30 // toroid -R 2 -r 1 -ns -3
31 // toroid -R 2 -r 1 -ns 3 -e 1
32 // toroid -R 2 -r 1 -ns 3 -e 1 -rs 1
33 // toroid -nphi 2 -ns 10 -e 1 -o 4
34 
35 #include "mfem.hpp"
36 #include <fstream>
37 #include <iostream>
38 
39 using namespace std;
40 using namespace mfem;
41 
42 static Element::Type el_type_ = Element::WEDGE;
43 static int order_ = 3;
44 static int nphi_ = 8;
45 static int ns_ = 0;
46 static double R_ = 1.0;
47 static double r_ = 0.2;
48 static double theta0_ = 0.0;
49 
50 void pts(int iphi, int t, double x[]);
51 void trans(const Vector &x, Vector &p);
52 
53 int main(int argc, char *argv[])
54 {
55  int ser_ref_levels = 0;
56  int el_type = 0;
57  bool dg_mesh = false;
58  bool visualization = true;
59 
60  OptionsParser args(argc, argv);
61  args.AddOption(&nphi_, "-nphi", "--num-elements-phi",
62  "Number of elements in phi-direction.");
63  args.AddOption(&ns_, "-ns", "--num-shifts",
64  "Number of shifts.");
65  args.AddOption(&order_, "-o", "--mesh-order",
66  "Order (polynomial degree) of the mesh elements.");
67  args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
68  "Number of times to refine the mesh uniformly in serial.");
69  args.AddOption(&R_, "-R", "--major-radius",
70  "Major radius of the torus.");
71  args.AddOption(&r_, "-r", "--minor-radius",
72  "Minor radius of the torus.");
73  args.AddOption(&theta0_, "-t0", "--initial-angle",
74  "Starting angle of the cross section (in degrees).");
75  args.AddOption(&el_type, "-e", "--element-type",
76  "Element type: 0 - Wedge, 1 - Hexahedron.");
77  args.AddOption(&dg_mesh, "-dm", "--discont-mesh", "-cm", "--cont-mesh",
78  "Use discontinuous or continuous space for the mesh nodes.");
79  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
80  "--no-visualization",
81  "Enable or disable GLVis visualization.");
82  args.Parse();
83  if (!args.Good())
84  {
85  args.PrintUsage(cout);
86  return 1;
87  }
88  args.PrintOptions(cout);
89 
90  // The output mesh could be hexahedra or prisms
91  el_type_ = (el_type == 0) ? Element::WEDGE : Element::HEXAHEDRON;
92  if (el_type_ != Element::WEDGE && el_type_ != Element::HEXAHEDRON)
93  {
94  cout << "Unsupported element type" << endl;
95  exit(1);
96  }
97 
98  // Determine the number of nodes in the cross section
99  int nnode = (el_type_ == Element::WEDGE)? 3:4;
100  int nshift = (ns_ >= 0) ? 0 : (nnode * (1 - ns_ / nnode));
101 
102  // Convert initial angle from degrees to radians
103  theta0_ *= M_PI / 180.0;
104 
105  // Define an empty mesh
106  Mesh *mesh;
107  mesh = new Mesh(3, nnode * (nphi_+1), nphi_);
108 
109  // Add vertices for a stack of elements
110  double c[3];
111  for (int i=0; i<=nphi_; i++)
112  {
113  c[0] = 0.0; c[1] = 0.0; c[2] = i;
114  mesh->AddVertex(c);
115 
116  c[0] = 1.0;
117  mesh->AddVertex(c);
118 
119  if (el_type_ == Element::HEXAHEDRON)
120  {
121  c[0] = 1.0; c[1] = 1.0;
122  mesh->AddVertex(c);
123  }
124 
125  c[0] = 0.0; c[1] = 1.0;
126  mesh->AddVertex(c);
127  }
128 
129  // Add Elements of the desired type
130  {
131  int v[8];
132  for (int i=0; i < nphi_; i++)
133  {
134  if (el_type_ == Element::WEDGE)
135  {
136  for (int j = 0; j < 6; j++) { v[j] = 3*i+j; }
137  mesh->AddWedge(v);
138  }
139  else
140  {
141  for (int j = 0; j < 8; j++) { v[j] = 4*i+j; }
142  mesh->AddHex(v);
143  }
144  }
145  }
146  mesh->FinalizeTopology();
147 
148  // Promote to high order mesh and transform into a torus shape
149  if (order_ > 1)
150  {
151  mesh->SetCurvature(order_, true, 3, Ordering::byVDIM);
152  }
153  mesh->Transform(trans);
154 
155  // Stitch the ends of the stack together
156  {
157  Array<int> v2v(mesh->GetNV());
158  for (int i = 0; i < v2v.Size() - nnode; i++)
159  {
160  v2v[i] = i;
161  }
162  // identify vertices at the extremes of the stack of prisms
163  for (int i=0; i<nnode; i++)
164  {
165  v2v[v2v.Size() - nnode + i] = (nshift + ns_ + i) % nnode;
166  }
167  // renumber elements
168  for (int i = 0; i < mesh->GetNE(); i++)
169  {
170  Element *el = mesh->GetElement(i);
171  int *v = el->GetVertices();
172  int nv = el->GetNVertices();
173  for (int j = 0; j < nv; j++)
174  {
175  v[j] = v2v[v[j]];
176  }
177  }
178  // renumber boundary elements
179  for (int i = 0; i < mesh->GetNBE(); i++)
180  {
181  Element *el = mesh->GetBdrElement(i);
182  int *v = el->GetVertices();
183  int nv = el->GetNVertices();
184  for (int j = 0; j < nv; j++)
185  {
186  v[j] = v2v[v[j]];
187  }
188  }
189  mesh->RemoveUnusedVertices();
190  mesh->RemoveInternalBoundaries();
191  }
192  if (order_ > 1)
193  {
194  mesh->SetCurvature(order_, dg_mesh, 3, Ordering::byVDIM);
195  }
196 
197  // Refine the mesh if desired
198  for (int lev = 0; lev < ser_ref_levels; lev++)
199  {
200  mesh->UniformRefinement();
201  }
202 
203  // Output the resulting mesh to a file
204  {
205  ostringstream oss;
206  if (el_type_ == Element::WEDGE)
207  {
208  oss << "toroid-wedge";
209  }
210  else
211  {
212  oss << "toroid-hex";
213  }
214  oss << "-o" << order_ << "-s" << ns_;
215  if (ser_ref_levels > 0)
216  {
217  oss << "-r" << ser_ref_levels;
218  }
219  oss << ".mesh";
220  ofstream ofs(oss.str().c_str());
221  ofs.precision(8);
222  mesh->Print(ofs);
223  ofs.close();
224  }
225 
226  // Output the resulting mesh to GLVis
227  if (visualization)
228  {
229  char vishost[] = "localhost";
230  int visport = 19916;
231  socketstream sol_sock(vishost, visport);
232  sol_sock.precision(8);
233  sol_sock << "mesh\n" << *mesh << flush;
234  }
235 
236  // Clean up and exit
237  delete mesh;
238  return 0;
239 }
240 
241 void trans(const Vector &x, Vector &p)
242 {
243  int nnode = (el_type_ == Element::WEDGE)? 3:4;
244 
245  double phi = 2.0 * M_PI * x[2] / nphi_;
246  double theta = theta0_ + phi * ns_ / nnode;
247 
248  double u = (1.5 * (x[0] + x[1]) - 1.0) * r_;
249  double v = sqrt(0.75) * (x[0] - x[1]) * r_;
250 
251  if (el_type_ == Element::WEDGE)
252  {
253  u = (1.5 * (x[0] + x[1]) - 1.0) * r_;
254  v = sqrt(0.75) * (x[0] - x[1]) * r_;
255  }
256  else
257  {
258  u = M_SQRT2 * (x[1] - 0.5) * r_;
259  v = M_SQRT2 * (x[0] - 0.5) * r_;
260  }
261 
262  p[0] = ( R_ + u * cos(theta) + v * sin(theta)) * cos(phi);
263  p[1] = ( R_ + u * cos(theta) + v * sin(theta)) * sin(phi);
264  p[2] = v * cos(theta) - u * sin(theta);
265 }
int visport
int main(int argc, char *argv[])
Definition: toroid.cpp:53
virtual void GetVertices(Array< int > &v) const =0
Returns element&#39;s vertices.
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:331
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
void Transform(void(*f)(const Vector &, Vector &))
Definition: mesh.cpp:12045
const Element * GetElement(int i) const
Return pointer to the i&#39;th element object.
Definition: mesh.hpp:1143
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:159
STL namespace.
void RemoveInternalBoundaries()
Definition: mesh.cpp:12203
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:1089
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition: mesh.hpp:1083
int AddVertex(double x, double y=0.0, double z=0.0)
Definition: mesh.cpp:1626
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
char vishost[]
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:10232
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Set the curvature of the mesh nodes using the given polynomial degree.
Definition: mesh.cpp:5635
Type
Constants for the classes derived from Element.
Definition: element.hpp:41
void trans(const Vector &x, Vector &p)
Definition: toroid.cpp:241
void RemoveUnusedVertices()
Remove unused vertices and rebuild mesh connectivity.
Definition: mesh.cpp:12096
int AddHex(int v1, int v2, int v3, int v4, int v5, int v6, int v7, int v8, int attr=1)
Definition: mesh.cpp:1757
void FinalizeTopology(bool generate_bdr=true)
Finalize the construction of the secondary topology (connectivity) data of a Mesh.
Definition: mesh.cpp:2945
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 GetNE() const
Returns number of elements.
Definition: mesh.hpp:1086
virtual int GetNVertices() const =0
RefCoord t[3]
int AddWedge(int v1, int v2, int v3, int v4, int v5, int v6, int attr=1)
Definition: mesh.cpp:1729
Vector data type.
Definition: vector.hpp:58
virtual void Print(std::ostream &os=mfem::out) const
Definition: mesh.hpp:2011
void pts(int iphi, int t, double x[])
double u(const Vector &xvec)
Definition: lor_mms.hpp:22
const Element * GetBdrElement(int i) const
Return pointer to the i&#39;th boundary element object.
Definition: mesh.hpp:1158
Abstract data type element.
Definition: element.hpp:28